Creation of this worksheet was supported by IT Academy program of Information Technology Foundation for Education (HITSA)

Exercises of Lab 8

The aim of the Lab is to derive an explicit finite difference method for solving an initial value problem of the heat equation in a bounded domain.

Let us consider the following problem: find $u$ such that \begin{gather}\frac{\partial u}{\partial t}(x,t)=\frac{1}{4}\frac{\partial^2 u}{\partial x^2}(x,t),\ x\in [-1,1],\ t\in (0,0.5]\\ u(-1,t)=1, u(1,t)=0,\ t\in (0,0.5]\\ u(x,0)=u_0(x),\ x\in [-1,1] \end{gather} where $u_0$ is a given function.

The procedure for deriving a finite difference approximation for the problem above consists of the following steps.

1. Choose a set of points at which we want to find approximate values of the unknown function. We define this set of points by dividing the interval $[-1,1]$ in $x$ direction into $n$ equal subintervals and the time interval $[0,0.5]$ into $m$ subintervals: we get points $(x_i,t_k)$, where $x_i=-1+i\frac{2}{n},\ i=0,\ldots,n$, $t_k=k\frac{0.5}{m}$. Since we know the values of the unknown function for $x=-1$, $x=1$ and for $t=0$, we have to determine approximate values $U_{ik}\approx u(x_i,t_k),\ i=1,\ldots,n-1,\ k=1,\ldots,m$, thus we have $m\cdot(n-1)$ unknowns (see the picture below for $n=3,m=4$).
2. In order to determine the values for $m\cdot (n-1)$ unknowns, we need $m\cdot (n-1)$ equations. We get those equations by writing down the differential equation at $m\cdot (n-1)$ points and then replacing the derivatives by approximations that use only the function values at points $(x_i,t_k),\ i=0,\ldots,n,\ k=0,\ldots,m$.
3. In order to get an explicit finite we should start using the equation at the points where we know the solution. For the current problem this means that we want to use points at the line, where the initial condition is given. So we use the equation at the points $(x_i,t_k),\ i=1,\ldots,n-1,\ k=0,\ldots,m-1$ and use the approximations \begin{align*} \frac{\partial u}{\partial t}(x_i,t_k)&\approx \frac{U_{i,k+1}-U_{ik}}{\Delta t}\ \text{(error}\leq const.\Delta t),\\ \frac{\partial^2 u}{\partial x^2}(x_i,t_k)&\approx \frac{U_{i-1,k}-2 U_{ik}+U_{i+1,k}}{\Delta x^2}\ \text{(error}\leq const.\Delta x^2). \end{align*} This gives us equations $$\frac{U_{i,k+1}-U_{ik}}{\Delta t}=\frac{1}{4}\frac{U_{i-1,k}-2 U_{ik}+U_{i+1,k}}{\Delta x^2}$$
4. Since in the equations corresponding to points on the line where we have data (initial condition, corresponding to $k=0$) the only unknown value is $U_{i,k+1}$, it makse sense to solve the equation for this value.

Exercise 1

Please find $a,b,c$ such that the we have $$U_{i,k+1}=a\,U_{i-1,k}+b\,U_{ik}+c\,U_{i+1,k}$$ and give a complete description of the method (how the matrix $U$ should be filled)

Check your answer (click to expand)

After following the steps of deriving the Explicit finite difference method for this problem, we get \begin{gather} U_{i0}=u_0(x_i),\ i=0,\ldots,n\\ U_{0k}=1,\ k=1,\ldots,m\\ U_{nk}=0,\ k=1,\ldots,m\\ U_{i,k+1}=\frac{1}{4}\frac{\Delta t}{\Delta x^2}\,U_{i-1,k}+(1-\frac{1}{4}\frac{2\Delta t}{\Delta x^2})\,U_{ik}+\frac{1}{4}\frac{\Delta t}{\Delta x^2}\,U_{i+1,k}, k=0,\ldots,m-1,\ i=1,\ldots,n-1 \end{gather}

Exercise 2

Write a function that for given values of $m$ and $n$ and for given function $u_0$ returns the values $U_{im},\ i=0,\ldots,n$ of the approximate solution obtained by explicit finite difference method. Test the correctness of your function in the case $m=100,n=10$ and $u_0(x)=\sin(\pi x)+\frac{1-x}{2}$, when the exact solution is $u(x,t)=e^{-\pi^2 t/4}\sin(\pi x)+\frac{1-x}{2}$. Hint: $\pi$ in python is np.pi

Solution

import numpy as np
def lab8solver(m,n,u0):
    xmin=??
    xmax=??
    T=???
    delta_t=??
    delta_x=? 
    #define values of x_i
    x=??
    #define matrix U with dimension (n+1)x(m+1)
    U=??
    #fill in the initial condition
    
    #use boundary conditions

    #the coefficients a,b,c in the formula for U[i,k+1]
    #do not depend on time, so compute before starting a time cycle
    a=
    b=
    c=
    #compute all other values 
    
    return U[:,m]
#define u0
def u0(x):
    return ??
m=100
n=10
approx_sol=lab8solver(m,n,u0)
print(approx_sol[:4])
#compare answers with values of exact solution

Check your answer (click to expand)

The first 4 values of the approximate solution should be [1. 0.72310846 0.51378348 0.41378348]
The first 4 values of the exact solution should be [1. 0.72882933 0.52304004 0.42304004]
The maximal absolute error of the solution should be 0.009256558574488039

Exercise 3

The total error caused by replacing exact derivatives with finite difference approximations is $O(\Delta t+\Delta x^2)$, which usually implies that the error of the approximate solution is of the same order. This means, that if we increase $m$ four times and $n$ two times, then the total error should be reduced approximately four times. Verify the convergence rate by computing the errors in the settings of the previous exercise for $m=4,16,64,256,1024$ and $n=2,4,8,16,32$.

#compute the errors and ratios of the consequtive errors

Check your answer (click to expand)

The errors are [3.56632987e-17 6.48611972e-02 1.53203711e-02 3.77115439e-03 9.39068948e-04]
The ratios of the errors are [5.49840277e-16 4.23365705e+00 4.06251496e+00 4.01584399e+00]

Exercise 4

It turns out that explicit methods may be unstable for certain choices of parameters $m$ and $n$. This means, that if $m$ and $n$ do not satisfy certain condition, the approximate solution may have arbitrarily large errors even when we let $m$ and $n$ to go to infinity. The sufficient condition of stability is that the coefficients $a$, $b$ and $c$ are all nonnegative. Repeat the computations of the previous exercise for $m=2,8,32,128$ and $n=10,20,40,80$ and compute the errors. Are the values of the approximate solution converging to the correct values? Find also the values of $a$, $b$ and $c$. Is the stability condition satisfied?

#comute the errors and ratios of errors

Check your answer (click to expand)

The errors are [1.22363261e-01 2.61927242e-02 1.88914222e+06 1.76749584e+75]
The ratios of the errors are [4.67165081e+00 1.38648768e-08 1.06882414e-69]

#compute a,b,c for each of the computation, print the values

Check your answer (click to expand)

The coefficients for the last computation are a = 1.5624999999999998 , b = -2.1249999999999996 , c = 1.5624999999999998