It has been a long time ago since my last Python tutorial, and I get a new computer since that, so.. first thing, installing again Python.
Ok, first problem. My tutorial for Python, Scipy and Matplolib is no longer valid. There is some problem with the versions and windows 8 and they don’t install properly. So I guess I need to do a new one.
I found that a good way of having Python and all the main scientific packages is to go for Anaconda.
Step 1. Go to the downloads web page. And look for the one to your system. In my case with Windows 8 and 64bits this is the one I choose:
Step 2. Just run the archive you downloaded. The installation is automatic. And voilà!
Step 3. Hello World. Once installed, when you run the “Launcher” this is what you see:
And for now, just launch spyder, the other features have no interest. It will look something like this:
wowow (I haven’t seen it before), It’s quite similar to Matlab, which I know very well. To start with it, type on the console:
Welcome to Spyder!
Now for the maths. Did you remember my previous post about how to solve the Rössler equation system? Check it back again, because we are going to do something similar.
Today we are going to focus on the diffusion equation.
Suppose we have some kind of medium (like water for example) and we deposit a drop of ink on it. The concentration of the ink at the point (x,y,z) at the instant t will be p(x,y,z,t). If we think that the ink diffuses by random movement of the atoms (water atoms and ink atoms), and there is no preferred direction (for example because of convection currents), then the evolution of the concentration of ink will be given by the diffusion equation:
The constant D determines how fast the diffusion takes place.
This equation is also called heat equation, because it also describes how heat diffuses on a solid body.
We are going to show how to use python to solve it numerically using the finite differences method.
First we need to discretize the domain in with we want to solve the function. So, instead of finding p(x,y,z,t) for every value of x, y, z and t, we are going to seek the solution only for the points:
determine how big is the separation is between the points where we are going to approximate the solution. Obviously, the smaller the factor the closer the points and the more accurate the solution (and more time to solve it).
Ok, so now that our domain divided into a grid, let’s try to write the diffusion equation using only points in the grid. According to the definition of derivative, it is true
To make it more readable we can use
If now instead of the real derivative we want the numerical approximation, we can take any of these:
There is no complication of going from total to partial derivatives in all these definitions (let’s do it from now on).
For the next step we are going to use a trick. Let’s take a point between two points of our grid
And write the central differences for two of those points:
Nice trick. Take your time to think about it and the grid.
Why to do this trick? Why not? We can select the grid we want. Doing so is a way of visualizing what is happening, instead we can start from the beginning with more points and nothing will change.
Using the trick, we can now write the second partial derivative as
And that will lead to
Now we can put together all the results and with the correct sub index. They look like…
OMG They look so complex! Not at all. If we introduce these expressions into the diffusion equation again and write the future value respect to the previous one…
And voilà! If we have the solution for all the points x, y, z at one instant of time, we can find the solution recursively for all the other instants!!!
Are we ready now for the Python code? No, not yet. Just two more things. we know that an initial condition is needed in order to solve the equation recursively, but what about boundary conditions?
There is mainly 3 types of boundary conditions.
Equal to a constant.
One example of this is that the boundaries of the problem are kept at a constante temperature (for instante, simulating an object inside a fridge or in an oven). This is achieved simple with
Typically, for the variable t instead of constants is more common using a density distribution only at the beginning.
This condition involves that nothing can cross the boundaries of the problem. This kind of condition can be for example ink in water inside a glass. The ink can difusse on the water, but cannot cross the walls of the glass.
To stablish this condition, we use the forward difference approximation
Making it equal to zero leads to:
And appliying the same trick to each boundary:
This boundary implies that the boundaries are connected, if we go enough time into the x direction we will eventually come back again to the starting position. This can be usefull when we are simulating a large periodic body and the same conditions must appear everywhere. It’s as simple as this:
Of course all this conditions can be mixed to create the kind of problem we need to solve.
Ok. Are we ready now? No, not yet, only one more thing.
When we deal with numerical approximations is crucial that the conditions we use meet the requirements to find a convergent solution. In this case, the size of the grid, the distant between approximations of the solutions, are the conditions we must look carefully. For the finite elements method with the approximations we have made, the requirement is:
Ok, so now we are ready.
The problem we are going to solve is a 2D square plate where 20 points inr andom positions have either +1 or -1 temperature, also randomly.
For the boundary conditions we are going to use free boundaries, no cosntant value, no periodic and no zero flux. This is just an example, is not correlated with any real problem I can imagine.
dt=0.002; Dz=1; dx=1; dy=1; #Those parameters fix the number of points in the grid. For isntance, the total x length of the plate will be dx*nx nx=120; ny=120; nt=50000; #We precharge the p matrix which will have inside the numerical solutions. p=zeros([nx,ny,nt]); #Here we set the boundary conditions for t (can be call initial conditions). It is going to be 20 hot or cold points in random positions. for f in range(20): p[round((nx-1)*rand()),round((ny-1)*rand()),0]=sign(rand()-0.51); for m in range(1,nt): #A simple implementation (but not quite eficient) will be iterating each point at a time. Because we have a matrix, we can operate with whole sections of the matrix at each time. #Basically, we take time slices and operate them as a whole. To use centered differences, we simply shift the matrix one element in the x direction or in the y direction. p[1:nx-1,1:ny-1,m]=p[1:nx-1,1:ny-1,m-1]+dt*Dz*((p[2:nx,1:ny-1,m-1]-2*p[1:nx-1,1:ny-1,m-1]+p[0:nx-2,1:ny-1,m-1])/np.power(dx,2)+(p[1:nx-1,2:ny,m-1]-2*p[1:nx-1,1:ny-1,m-1]+p[1:nx-1,0:ny-2,m-1])/np.power(dy,2)); #Finally we plot several iterations fig = figure() for g in range(18): subplot(3,6,g+1) axis([0,nx,0,ny]) pcolor(p[:,:,round(g*(nt-1)/17)],vmin=-0.001,vmax=0.001,cmap='jet') axis('off')
It took more time to explain the method than the code itself. Amazing.
This is the output from the code.
Nice. Pretty nice.
And that was all, hope you like it.
Ok, here it is an animation of the evolution of the heat for that particular solution.