Partial Differential Equations in Python.
Diffusion Equation in Python.
What are Partial Differential Equations (PDEs)?
Partial Differential Equations are the equations that have more than one independent variables, an unknown function which depends on those variables and partial derivatives of the unknown function with respect to the independent variables. Please, read the references for more.
The general for of a PDE is:
By considering only the first three coefficients, A B and C, we can determine what equation are we dealing with or what problem are we solving.
- If 𝐵²−4𝐴𝐶 < 0, we have an elliptic PDE. On this king of problems, the Laplace’s Equation us used.
The two derivatives of this equations are the derivatives of space x² and y², no time derivative.
- If 𝐵²−4𝐴𝐶 >0, then we have a hyperbolic PDE, where the Wave Equation is used
The two derivatives of this equation are the Time in second order t² and a space derivative in second order y².
- If 𝐵²−4𝐴𝐶 = 0, then we have a parabolic PDE, and the Diffusion Equation is used.
Just like the hyperbolic PDE, the Diffusion looks similar, but the time here is in first order and space is second order.
Computers and Differentials
- Forward Difference
Although computers are very good at math, they do not understand differential equations. In order to tell the computer to solve differential equations, we need to discretize the equation. Consider the following:
U=u(x, y)
The partial expansion in terms of x is
By discarding the second order terms, we will end up with a very approximate formula:
This is a famous Euler Method scheme, which has been seen in ODE (Ordinary Differential Equations). It is called Forward Difference in some finite difference.
Take a look at the following example:
The forward difference is:
Backward Difference
By looking at the forward difference formula, we notice that the the forward difference is found by subtracting next step uᵢ₊₁,ⱼ with the current step uᵢ,ⱼ. Intuitively, we can think of backward difference.
The backward difference can be found by subtractin current step uᵢ,ⱼ with the previous step uᵢ-1,ⱼ
Central Difference
Consider the following
Discarding the third order onward, we will end with an approximation
Similarly to forward and backward difference , but this time we “stand” at the center, meaning that on the left side 𝑢(𝑥+Δ𝑥) we look forward and the right side 𝑢(𝑥−Δ𝑥) we look backward. This is the Central Difference because you stand at the center and look forward differences and backwards as well/.
In terms of indexes, we will have the combined formula.
Diffusion Equation
To solve the Diffusion Equation in 1D in a finite system u(x, t), where t is time, we write:
… for 0 ≤ 𝑥 ≤ 𝐿, and D is the thermal diffusivity D = √(k/Cp), where:
- K: thermal conductivity
- p (greek ro): density of the material
- C: the specific heat
Looking at X, we notice that it starts from 0 and ends at L, so 0 and L are the boundaries of X, where L is the length of a rod.
To solve such equations, an initial condition is needed, e.g 𝑢(𝑥,0) = 𝑢₀(𝑥), meaning that the temperature at the time 0 is 𝑢₀(𝑥). The initial condition is always some function, which represent the temperature of all points at time 0, later we see how it evolves.
Two boundary conditions are needed as well for solving the equation, where the boundaries will be fixed over time. 𝑢(0,𝑡)=𝑢₁(𝑡), which is the left boundary and 𝑢(𝐿,𝑡)=𝑢₂(𝑡) for the right boundary.
The boundary conditions are functions of time, which means they will be constant over time.
Note: This is a very simple form of Diffusion Equation because radiation, cooling and other factors are not considered. Or the rod 0≤𝑥≤𝐿 is isolated from the environment.
Writing the equation in discretized form (forward difference and central difference):
…where the total error is O(k + h²)
Note: In the first formula we took j+1 because we assume j is the time index, bun in second formula i+1 because i is space in second order.
The discretized formula will look like this..
Simplifying the formula, we take 𝐷=√(k𝐶/𝑝)=1 for simplicity. Since D = 1, we take the left denominator k and we factor by right denominator ℎ² which gives us the ratio 𝑟=𝑘/ℎ², and the simplified form will look look like below:
Coding Diffusion in 1D
Now that we discretize the formula, we leave to the computer to do the math.
Coding summary.
First import numpy and matplotlib. Create a numpy linspace array for the X values and a np.zeros with the exact values as X for the U function. Declare the initial condition and boundary conditions. Lastly, loop for the time, say 100, and then for the space starting from 1 to L-1 because we assign the boundary conditions and they are fixed.
The above code will produce the following plot.
What about Diffusion in 2D?
In 1D we saw that ∂u\∂t = ∂²𝑢\∂𝑦², now the same formula applies but we add another second order space derivative.
Moving on, descretize the formula:
…where:
- i: x-index
- j: y-index
- t: time index
Simplifying the formula, we end up with with:
Coding Diffusion in 2D
- Coding summary
Import numpy and matplotlib as we did in 1D. Also import seaborn for a bit style. For plotting the results, we will use contourf function in matplotlib, so x and y arrays will not be needed. Define u as np.zeros with as shape of 100x100. Define an initial condition but since the boundaries will be 0 and we have zeros grid, we will not define the boundaries.
To make the figure more fun, lets do it with an animation in 2D first and then 3D.
2D animation code
Results as below…
3D animation
- Code summary
As the 2D, the 3D will be very similar. But in order to create a 3D shape with countourf3D in matplotlib, we actually need x and y from a np.meshgrid. In contourf3D we will pass X, Y and u.
Results as bellow…
Conclusion
PDEs are equations that have many independent variables. They have an unknown function which depends on the independent variables and partial derivatives of the unknown function. At its general form, we can determine what problem are we solving. If 𝐵²−4𝐴𝐶<0, we have an elliptic PDE. If 𝐵²−4𝐴𝐶>0, we have a hyperbolic PDE. And if 𝐵²−4𝐴𝐶=0 then we have a parabolic PDE. Computer are way better than humans when it comes to solving arithmetic problems, but they don’t understand differential equations. So in order to tell to the computer how to solve them, we must first discretize the equation.
In general, PDEs are just mathematical tools to model engineering systems such as heat transfer, as we did in the examples above, hydraulic flow, electrical circuits and more…
