1.9. Treatment of Partial Differential Equations: a start

Partial Differential Equations (PDEs) are characterized by having more than one independent variable. They can be classified as the ODEs, in linearity or non linearity, first order, second order or higher. See the following widely used PDEs:

\[ u = -k_x \frac{\partial \phi}{\partial x}; \text{ } v = -k_y \frac{\partial \phi}{\partial y}; \text{ Darcy's law -> PDE system of first order} \]

\[ \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = 2xy \text{ Poisson's equation -> PDE second order} \]

\[ \frac{du}{dt} = \nu \frac{d^2u}{dx^2} \text{ diffusion equation -> PDE second order} \]

\[ \frac{\partial \phi}{\partial t} + c \frac{\partial \phi}{\partial x} = 0 \text{ Convection -> PDE first order} \]

As you can see, there are PDEs that depend only on space variables (e.g. poisson equation) and those that depend on time and space (e.g. diffusion equation). We will now mainly treat the latter using the method of lines. Let’s look at a new notation. Consider the variable \(u\) at the location \(x=i\) and time \(t=j\). It can be written as:

\[ u(x=i,t=j)=u^j_i \]

Method of Lines

This method transforms the PDE into a system of ordinary differential equations by replacing the spatial derivatives with a numerical approximation. Let’s see this transformation in the diffusion equation using the dependent variable temperature \(T\) and the independent variable \(x\) for space:

\[ \frac{dT}{dt} = \nu \frac{d^2T}{dx^2} \]

Clearly \(T(x,t)\) varies in time and space. Due to the first derivative, an initial condition is needed: the initial temperature defined in the whole space of interest. Due to the second derivative, two boundary conditions are needed. The spatial second derivative can be approximated with any method, in this example a central difference is used:

\[ \left. \frac{dT}{dt} \right|_{x_i,t} = \frac{\nu}{\Delta x ^2} \left( T^t_{i-1} - 2T^t_{i} + T^t_{i+1} \right ) \text{ for } i=1,2,...n-2,n-1 \]

The counter of that equation does not start at i=0 and it does not end at \(n\) due to the BC. For constant Dirichlet conditions, \(T^t_{0/n}=T_{left/right}\) at all times. Let’s consider a problem with a grid of 5 points with only three unknowns (again due to the BC), the system of ODEs becomes:

\[ \frac{dT^t_1}{dt} = \frac{\nu}{\Delta x ^2} \left( T_{0} - 2T^t_{1} + T^t_{2} \right ) \]

\[ \frac{dT^t_2}{dt} = \frac{\nu}{\Delta x ^2} \left( T^t_{1} - 2T^t_{2} + T^t_{3} \right ) \]

\[ \frac{dT^t_3}{dt} = \frac{\nu}{\Delta x ^2} \left( T^t_{2} - 2T^t_{3} + T_{4} \right ) \]

Solution with Forward Euler for the time discretization

Then, the time derivative is approximated, with any method, here we use the explicit method Forward Euler:

\[ \frac{T^{j+1}_i-T^{j}_i}{\Delta t} = \frac{\nu}{\Delta x ^2} \left( T^j_{i-1} - 2T^j_{i} + T^j_{i+1} \right ) \text{ for } i=1,2,...n-2,n-1; \text{ for } j=0,1,2,...,m \]

To approximate the temperature for \(t=1\), the system of equations starts at \(j=0\), as it is an explicit method, we can already solve for the next time step:

\[ T^1_1 = T^0_1 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{0} - 2T^0_{1} + T^0_{2} \right ) \]

\[ T^1_2 = T^0_2 + \frac{\nu \Delta t}{\Delta x ^2} \left( T^0_{1} - 2T^0_{2} + T^0_{3} \right ) \]

\[ T^1_3 = T^0_3 + \frac{\nu \Delta t}{\Delta x ^2} \left( T^0_{2} - 2T^0_{3} + T_{4} \right ) \]

The solution for \(T^1_{1,2,3}\) is directly found because \(T^0_{1,2,3}\) are known thanks to the initial condition.

Solution with Backward Euler for the time discretization

If instead of an explicit Forward Euler, we use an implicit scheme like Backward Euler, the system of ODEs becomes:

\[ T^1_1 = T^0_1 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{0} - 2T^1_{1} + T^1_{2} \right ) \]

\[ T^1_2 = T^0_2 + \frac{\nu \Delta t}{\Delta x ^2} \left( T^1_{1} - 2T^1_{2} + T^1_{3} \right ) \]

\[ T^1_3 = T^0_3 + \frac{\nu \Delta t}{\Delta x ^2} \left( T^1_{2} - 2T^1_{3} + T_{4} \right ) \]

Separating the unknowns:

\[ T^1_1 - \frac{\nu \Delta t}{\Delta x ^2} \left( - 2T^1_{1} + T^1_{2} \right ) = T^0_1 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{0} \right ) \]

\[ T^1_2 - \frac{\nu \Delta t}{\Delta x ^2} \left( T^1_{1} - 2T^1_{2} + T^1_{3} \right ) = T^0_2 \]

\[ T^1_3 - \frac{\nu \Delta t}{\Delta x ^2} \left( T^1_{2} - 2T^1_{3} \right )= T^0_3 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{4} \right ) \]

This can be written as \(Ay=b\) and can be solved in a similar way as the second order ODE of the previous chapter:

\[\begin{split} \begin{bmatrix} (1 + 2 \frac{\nu \Delta t}{\Delta x ^2}) & - \frac{\nu \Delta t}{\Delta x ^2} & 0 \\ - \frac{\nu \Delta t}{\Delta x ^2} & (1 + 2 \frac{\nu \Delta t}{\Delta x ^2}) & - \frac{\nu \Delta t}{\Delta x ^2} \\ 0 & - \frac{\nu \Delta t}{\Delta x ^2} & (1 + 2 \frac{\nu \Delta t}{\Delta x ^2}) \\ \end{bmatrix} \begin{bmatrix} T^1_1 \\ T^1_2 \\ T^1_3 \\ \end{bmatrix} = \begin{bmatrix} T^0_1 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{0} \right )\\ T^0_2 \\ T^0_3 + \frac{\nu \Delta t}{\Delta x ^2} \left( T_{4} \right ) \\ \end{bmatrix} \end{split}\]

For the \(T*2_{1,2,3}\) the solution can be found in the same way given that now \(T*1_{1,2,3}\) are known.

Procedure Summary

Some of the following steps can be re ordered without affecting the solution. The following structure is proposed:

1. Discretization of the spatial derivative with a numerical approximation

2. Discretization of the time derivative with a numerical approximation

3. Parameter definition

4. Grid creation

5. Define Initial conditions

6. Define Boundary conditions

7. Building a system of equations according to the discretization

8. Solve the system

9. Go back to 7 until the final time step, \(m\), is reached.