Lessons learned about FEM with SciPy

I am just sitting in a train doing a programming exercise about the Finite Element Method.

Details about the examples can be found on the homepage of Prof. Sormann at http://itp.tugraz.at/LV/sormann/AKNumPhysik/

The first problem is solving the stationary heat equation with FEM.

Now how do we plot this? At first I tried using matplotlib which was not satisfactory because it cannot do color gradients across triangles in a mesh grid. Then I discovered the wonderful MayaVi package where you can do it as simple as


from enthought.mayavi import mlab
mlab.triangular_mesh(x, y, z, triangle_indices, temperature)
mlab.show()


That's all! To add fancy stuff like a wireframe mesh, a colorbar and axes I did something like that before the final mlab.show()

mlab.triangular_mesh(x, y, z, triangle_indices, color=(0,0,0),

                     line_width=1.,representation='wireframe')

mlab.view(0., 0.)

mlab.move(0.,-.5,0.)

mlab.xlabel("x")

mlab.ylabel("y")

mlab.axes(extent=[0., 4., 0., 4., 0., 0.], nb_labels=5)

mlab.colorbar(title="Temperatur", orientation='vertical')

And here is the result:

Nice, huh?


If I have time, I'll keep you updated about the algorithm, the progress for the time-dependent heat conduction and more!

Equivalent ODE integrators in Matlab and SciPy

Here is a short summary of my investigations on how ODE integrators from SciPy compare with the ones provided by Matlab:

• SciPy provides several ODE solvers based on Fortran routines from ODEPACK
• The default SciPy solver is LSODA and provided via scipy.integrate.odeint(f)
• There are equivalent solvers for Matlab's ode45 and ode15s.
• Use scipy.integrate.ode(f).set_integrator(...)
• ode45: 'dopri5'
• ode15s: 'vode', method='bdf', order=15

Something else of interest : NodePy looks like a tool to experiment with various ODE solving algorithms and offers more flexibility than either Matlab or SciPy.

In conclusion one can say that if you know what you are doing, there are plenty of good integrators available for Python and it's easy to plug in your own Fortran or C code.

Links:
[1] http://www.scipy.org/NumPy_for_Matlab_Users
[2] http://www.mathworks.com/help/techdoc/ref/ode23.html
[3] https://github.com/scipy/scipy/blob/master/scipy/integrate/ode.py
[4] http://web.kaust.edu.sa/faculty/davidketcheson/NodePy/index.html

[5] https://computation.llnl.gov/casc/odepack/odepack_home.html

MathJax in Blogger

Hurray - now I can finally use $\LaTeX$ in Blogger, thanks to a post on http://mnnttl.blogspot.com/ !
Look at the Schrödinger Equation for example:

$$
\mathrm{i}\hbar\frac{\partial}{\partial t} |\,\psi (t) \rangle = \hat{H} |\,\psi (t) \rangle
$$

Now that's nice... but what about the Lagrange density? ... wait, here it is:

$$
\mathcal{L}\left(\psi, \mathbf{\nabla}\psi, \dot{\psi}\right) = \mathrm i\hbar\, \frac{1}{2} (\psi^{*}\dot{\psi}-\dot{\psi^{*}}\psi) - \frac{\hbar^2}{2m} \mathbf{\nabla}\psi^{*} \mathbf{\nabla}\psi - V( \mathbf{r},t)\,\psi^{*}\psi
$$

And now for some more complicated stuff:

$$
\begin{align}
\mathcal{L}_\mathrm{QCD}(q, A)
& = \bar{q}\left(i \gamma^\mu D_\mu - m \right) q - \frac{1}{4}F^a_{\mu \nu} F^{\mu \nu}_a \\
& = \bar{q} (i \gamma^\mu \partial_\mu - m) q + g \bar{q} \gamma^\mu T_a q A^a_\mu - \frac{1}{4}F^a_{\mu \nu} F^{\mu \nu}_a \\
\end{align}
$$

We should definitely use that to teach physics!

Blog

Hm... heute dachte ich mir, ich brauche in Blog... wozu auch immer :-)