Algorithm suggestion for PDE - example: heat equation

Question

I want to solve the PDE equation numerically. For this, I started my study with something simple; heat equation $$ \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial^2 x} $$ with the initial condition $$ u(x,0)=1\qquad(-1<x<1) $$ I assumed the constant is set to zero.

My next step was to turn this into discrete equation. I used $$ \frac{\partial u}{\partial t}=\frac{u(t+\Delta t)-u(t)}{\Delta t},\qquad \frac{\partial^2 u}{\partial x^2}=\frac{u(x+2\Delta x)+u(x)-2u(x+\Delta x)}{(\Delta x)^2} $$

Arranging all these, $$ \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t} = \frac{u(x+2\Delta x,t)+u(x,t)-2u(x+\Delta x,t)}{(\Delta x)^2} $$ $$ u(x,t+\Delta t)=u(x,t)+\frac{\Delta t}{(\Delta x)^2}\left( u(x+2\Delta x,t)+u(x,t)-2u(x+\Delta x,t) \right) $$ so, at arbitrary $x$, after time $\Delta t$, the value will be function of $u$ at $x+2\Delta x$ and $x+\Delta x$, before $\Delta t$.

This is my algorithm, and I wrote a code for Matlab

clc; clear all;
pt=1000;                    %%% number of points
xsta=-5;                    %%% x start
xend=5;                     %%% x end
x=linspace(xsta,xend,pt);   %%% creating x space -5 ~ 5
t=linspace(0,2,1000);       %%% creating t space 0 ~ 2 sec

dx=abs(x(2)-x(1))/pt;       %%% x interval
dt=abs(t(2)-t(1))/pt;       %%% t interval

u=zeros(pt,pt);             %%% creating u function

u0=zeros(1,pt);             %%% initial value
for i=1:pt
    if x(i)>=-1 && x(i)<=1
        u0(i)=1;
    end
end                         %%% initial value

u(1,:)=u0;                  %%% inserting initial value to t=0;

for i=2:pt                  %%% t step. it starts from 0+dt
    for j=1:pt-2;           %%% x step. it ends at xend - dx*2 due to algorithm
        u(i,j)=u(i-1,j)+((  u(i-1,j+2) + u(i-1,j) -2*u(i-1,j+1)  )/dx)*dt;
    end
end

figure();
plot(x,u(150,:))            %%% plotting at t=0.2983 sec

But when I plot, It gives weird function, and does not look like heat equation.

Moreover, the scale is off by $10^{36}$

I am writing my code for first time, so I don't know what is wrong with my code.

Do you have any suggestion of better algorithm for this PDE?

Also, this is just start. I am going to use more complicated, coupled PDE for my project. Is there any good algorithm for this?

Answer 1

While i agree with Wolfgang that its best to choose Crank-Nicolson time stepping, i think it is incorrect to assert explicit time stepping results in 'incredibly small' time steps or even to suggest your script is not functioning because of explicit time stepping. You can make it work just fine, but you need to get your discretization correct and to implement the boundary conditions properly.

Problems:

First, your discretization of your spatial second derivative is incorrect. I like to do this using the Finite Volume method because you can not go wrong. You define an average temperature over a control volume $[x-\frac{1}{2}\Delta x,x+\frac{1}{2}\Delta x]$ where $x$ is the node at the center of the volume element with size $\Delta x$: $$\bar{u} = \frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} u dx$$

Now you integrate the heat equation: $$\partial_{t}u=a\partial_{x}^{2}u$$ over the control volume: $$\frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} \partial_{t}u dx=a\frac{1}{\Delta x}\int_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x} \partial_{x}^{2}u dx$$ $$\partial_{t}\bar{u} =\frac{a}{\Delta x}\left.\partial_{x}u \right|_{x-\frac{1}{2}\Delta x}^{x+\frac{1}{2}\Delta x}$$

Edit: The result here is that the temperature in the control volume is spatially uniform, i.e. the average temperature $\bar{u}$ at node $x$. For $\Delta x$ small enough, this become increasingly more accurate. The change in time of $\bar{u}$ then becomes dependent of the incoming and outgoing flux of heat at the boundaries located at $x-\frac{1}{2}\Delta x$ and $x+\frac{1}{2}\Delta x$. The choice of these coordinates make this a staggered grid where the nodes do not align with the boundaries.

using the second-order accurate discretization of the derivate, we get: $$\partial_{t}\bar{u} = \frac{a}{\Delta x}\left[\frac{\bar{u}\left(x+\Delta x\right)-\bar{u}\left(x\right)}{\Delta x}-\frac{\bar{u}\left(x\right)-\bar{u}\left(x-\Delta x\right)}{\Delta x}\right] = \frac{a}{\Delta x^{2}}\left[\bar{u}\left(x+\Delta x\right)-2\bar{u}\left(x\right)+\bar{u}\left(x-\Delta x\right)\right] $$

which is different from your discretization. Edit: since the averaged temperature at node $x$ is used exclusively in these equations i will drop the 'bar' from $\bar{u}$

Edit: For abritrary coordinate $\tilde{x}$, second-order accuracy is shown by Taylor expansion: $$u\left(\tilde{x}-\frac{\Delta x}{2}\right)=u\left(\tilde{x}\right)-\frac{\Delta x}{2}\left.\frac{du}{dx}\right|_{\tilde{x}}+\frac{1}{2}\left(\frac{\Delta x}{2}\right)^{2}\left.\frac{d^{2}u}{dx^{2}}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$

$$u\left(\tilde{x}+\frac{\Delta x}{2}\right)=u\left(\tilde{x}\right)+\frac{\Delta x}{2}\left.\frac{du}{dx}\right|_{\tilde{x}}+\frac{1}{2}\left(\frac{\Delta x}{2}\right)^{2}\left.\frac{d^{2}u}{dx^{2}}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$

Subtracting to cancel the second derivative gives: $$u\left(\tilde{x}+\frac{\Delta x}{2}\right)-u\left(\tilde{x}-\frac{\Delta x}{2}\right)=\Delta x\left.\frac{du}{dx}\right|_{\tilde{x}}+O\left(\Delta x\right)^{3}$$

rewriting to an expression for the first derivative: $$\left.\frac{du}{dx}\right|_{\tilde{x}}=\frac{u\left(\tilde{x}+\frac{\Delta x}{2}\right)-u\left(\tilde{x}-\frac{\Delta x}{2}\right)}{\Delta x}+O\left(\Delta x\right)^{2}$$

substituting in the averaged heat equation and filling in the $\tilde{x}=x-\frac{1}{2}\Delta x$ and $\tilde{x}=x+\frac{1}{2}\Delta x$ respectively gives the discretized equation found before.

Secondly, you need to define boundary conditions. Since you didn't specify i am going to assume you assume zero-gradient boundary conditions. Since $x$ is at the center of the node, you can view it as a staggered grid. This means that the boundaries of the domain are in between nodes at $x_{b0}=x_1-\frac{1}{2}\Delta x$ and $x_{bf}=x_N+\frac{1}{2}\Delta x$, where $x_{1}$ and $x_{N}$ are the first and $N$th node respectively. We can set the gradients at $x_{b0}$ and $x_{bf}$ by introducing 'ghost' nodes outside the domain at $x_0$ and $x_{N+1}$ and then the gradients are defined as: $$\partial_{x}u\left(x_{b0}\right)=\frac{u\left(x_{1}\right)-u\left(x_{0}\right)}{\Delta x}=0\quad\partial_{x}u\left(x_{bf}\right)=\frac{u\left(x_{N+1}\right)-u\left(x_{N}\right)}{\Delta x}=0$$

So simply setting $u\left(x_{0}\right)=u\left(x_{1}\right)$ and $u\left(x_{N+1}\right)=u\left(x_{N}\right)$ will ensure zero-gradients at the boundaries.

Now we are in a position to solve the equations. Edit: For each node we solve the discretized heat equation and use the information at the boundaries for the first and last node. For simplicity each node is identified by a subscript $i$, i.e. $u_i=u\left(x_i\right)$:

$$\partial_{t}\bar{u}_{i}=\frac{a}{\Delta x^{2}}\left[\bar{u}_{i+1}-2\bar{u}_{i}+\bar{u}_{i-1}\right]$$

written out explicitly:

$$\partial_{t}\bar{u}_{1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{2}-2\bar{u}_{1}+\bar{u}_{0}\right]$$ $$\partial_{t}\bar{u}_{2} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{3}-2\bar{u}_{2}+\bar{u}_{1}\right]$$ $$\vdots $$ $$\partial_{t}\bar{u}_{N-1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{N}-2\bar{u}_{N-1}+\bar{u}_{N}\right]$$ $$\partial_{t}\bar{u}_{N} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{N+1}-2\bar{u}_{N}+\bar{u}_{N-1}\right]$$

Edit: Here $u_0$ and $u_{N+1}$ are the values of the temperature at the ghost nodes previously determined from the boundary conditions. We substitute those in to get:

$$\partial_{t}\bar{u}_{1} = \frac{a}{\Delta x^{2}}\left[\bar{u}_{2}-\bar{u}_{1}\right]$$ $$\partial_{t}\bar{u}_{N} = \frac{a}{\Delta x^{2}}\left[-\bar{u}_{N}+\bar{u}_{N-1}\right]$$

Next, we define a system of equations with vector: $$\vec{u}=\begin{bmatrix}u\left(x_{1}\right) & u\left(x_{2}\right) & \cdots & u\left(x_{N-1}\right) & u\left(x_{N}\right)\end{bmatrix}^{T}$$

which allows us to write these equations as a system of equations: $$\partial_{t}\begin{bmatrix}\bar{u}_{1}\\ \bar{u}_{2}\\ \vdots\\ \bar{u}_{N-1}\\ \bar{u}_{N} \end{bmatrix}=\frac{a}{\Delta x^{2}}\begin{bmatrix}-1 & 1\\ 1 & -2 & 1\\ & \ddots & \ddots & \ddots\\ & & 1 & -2 & 1\\ & & & 1 & -1 \end{bmatrix}\begin{bmatrix}\bar{u}_{1}\\ \bar{u}_{2}\\ \vdots\\ \bar{u}_{N-1}\\ \bar{u}_{N} \end{bmatrix}$$

so we get the linear system: $$\partial_t\vec{u}=\boldsymbol{A}\cdot\vec{u}$$

where: $$A=\frac{a}{\Delta x^{2}}\left[\begin{array}{ccccc} -1 & 1\\ 1 & -2 & 1\\ & \ddots & \ddots & \ddots\\ & & 1 & -2 & 1\\ & & & 1 & -1 \end{array}\right]$$

All rows are described by our discretized equation, except in the first and last column where the values for the 'ghost' nodes outside of the domain are used to account for the boundary conditions.

Let's not worry about discretization of the time derivative, we will simply let matlab integrate this system using an appropriate 'ode' function. This has as advantage that matlab choses an proper timestep for discretization which makes sure that the integration is stable.

Below you will find an implementation of the above linear system and the result.

Code and result:

function main()
    clc, clear all, close all;

    a = 1;                     %% thermal diffusivity
    pt=100;                    %%% number of points
    xsta=-5;                    %%% x start
    xend=5;                     %%% x end
    x = linspace(xsta,xend,pt);   %%% creating x space -5 ~ 5

    % create discretization matrix
    A = zeros(pt,pt);
    A(1,1) = -1; 
    A(1,2) = 1;
    for i = 2:pt-1
        A(i,i-1) = 1;
        A(i,i) = -2;
        A(i,i+1) = 1;
    end
    A(pt,pt-1) = 1; 
    A(pt,pt) = -1;

    % initial condition
    u0 = zeros(pt,1);
    u0(abs(x)<1) = 1;

    function dudt = pde(t, u, Fo_x)
        dudt = Fo_x*A*u;
    end

    nt = 5; % increase nt if you want more time points (not needed for accuracy)
    Fo = 1; % Fo is a dimensionless diffusion coefficient, a*t/l^2
    Fo_x = Fo*(pt-1); % Fo_x is diffusion coefficient based on grid
    tspan = linspace(0,1,nt);
    [dimt,u] = ode45(@pde, tspan, u0, [], Fo_x);

    figure, plot(x,u,'o')
    legendCell = cellstr(num2str((dimt*(xend-xsta)^2/a), 't=%-d s'));
    legend(legendCell)
    xlabel('spatial coordinate, x')
    ylabel('temperature, u')
end

Answer 2

You really don't want to solve the heat equation with an explicit time stepping scheme. You need to choose the time step so incredibly small that you won't make any progress towards the end time.

(Explanation in lecture 27: http://www.math.tamu.edu/~bangerth/videos.html .)

Answer 3

So you can solve this problem with an Explicit time stepping scheme, as nluigi mentioned, but it comes at a cost just as Wolfgang stated.

In the problem statement, the problem is formulated using an Explicit Euler time stepping scheme with (arguably) a second order central difference.

The stability region of Explicit Euler can be found to be:

$$\Delta t \leq \frac{2}{|\lambda|_{max}}$$

Where $|\lambda|_{max}$ is the largest modulus of the eigenvalues taken from the differential operator $A$ that allows the following discretization to be true:

$$ \frac{\partial \textbf{u} }{\partial t} = A \textbf{u}$$

Using the Second Order Central Difference scheme for space, and assuming Dirichlet boundary conditions, one can find the following expression for the largest eigenvalue modulus:

$$ |\lambda|_{max} = \frac{2 C}{\Delta x^2}$$

where $C \gt 1$ and $\Delta x$ represents the distance between points in the spacial discretization.

When you substitute this result into the stability result for Explicit Euler, you find:

$$ \Delta t \leq \frac{\Delta x^2}{C}$$

This result implies that as you shrink the distance between points, you must shrink the time step an order of magnitude more to ensure the solution remains stable. This is often unrealistic for higher precision simulations. This is the problem that Wolfgang was mentioning and why you should consider another scheme. But as long as you stay within the stability region, you can integrate with an Explicit scheme just fine, as nluigi mentioned.

To work around this, it's common to either use an implicit time stepping scheme or use an adaptive explicit scheme (which is what ode45 in Matlab uses).