We considered the model problem \[ u_{xx} + u_{yy} + u^2 = f(x,y) \] on \[ 0 \le x \le 1, \,\, 0 \le y \le 1 \] subject to \[ u(0,y) = u(1,y) = u(x,0) = u(x,1) = 0 \]
We finite difference the system using \( O(h^2) \) finite differences as previously. To implement SOR, we appeal to the general relation from the notes (for Gauss-Seidel relaxation): \[ u_{i,j}^{(n+1)} = u_{i,j}^{(n)} - r_{i,j}^{(n)} \left[ \frac{\partial F^h_{i,j}}{\partial u_{i,j}} \Bigg\vert_{u_{i,j}=u_{i,j}^{(n)}}\right]^{-1} \] where \( F^{h}_{i,j} = 0 \) is the difference equation for point \( (i,j) \) on the mesh, and \( r_{i,j}^{(n)} \) is the residual with respect to that equation. Here we have \[ \frac{\partial F^h_{i,j}}{\partial u_{i,j}} \Bigg\vert_{u_{i,j}=u_{i,j}^{(n)}} = -4h^{-2} + 2 u_{i,j} \]
The corresponding part of function sor2d_nl.m is
% Update unknown and accumulate residual ...
rij = (u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - ...
4*u(i,j)) / h^2 + u(i,j)^2 - f(i,j);
rnorm = rnorm + rij^2;
% Compute GS value ...
uij_gs = u(i,j) - rij / (-4.0 / h^2 + 2.0 * u(i,j));
% Overrelaxation step ...
u(i,j) = omega * uij_gs + (1.0 - omega) * u(i,j);
Given two difference solutions \( u^{h} \) and \( u^{2h} \) and the expansion
\[ \lim_{h\to0} u^h(x,y) = u(x,y) + h^2 e_2(x,y) +O(h^4) \] we seek a linear combination of the solutions that eliminates the \( O(h^2) \) error term. We have
\[ \alpha u^{h} + \beta u^{2h} = u(x,y) + O(h^4) \] from which we get \[ \alpha + \beta = 1 \] \[ \alpha + 4 \beta = 0 \] Solving, we find \[ \alpha = \frac 43 \quad\quad \beta = -\frac 13 \] so our extrapolation formula is \[ \frac 43 u^h - \frac 13 u^{2h} = u + O(h^4) \]
When I perform extrapolation of levels 6 and 7 I find the following
Level 6 error: 0.000106666
Level 7 error: 2.68702e-05
Extrapolated error: 3.86917e-09
so for this equation and particular levels of discretization, extrapolation reduces the error by almost four orders of magnitude.
Recall that we are solving
\[ u_t(x,t) = u_{xx} \] (unit diffusion constant) on \[ 0 \le x \le 1 \quad\quad t \ge 0 \] subject to \[ u(x,0) = u_0(x) \] and \[ u(0,t) = u(1,t) = 0 \] The scripts that you will use below are configured so that the initial data is \[ u_0(x) = \sin(2\pi x) \] corresponding to a solution \[ u(x,t) = e^{-4\pi^2t}\sin(2\pi x) \]
This scheme has \( O(\Delta t, \Delta x^2) \) truncation error.
DOWNLOAD SOURCE
EXERCISE
Run the script tdiff_1d_ftcs.m. What do you observe? Now modify
tdiff_1d_ftcs.m so that the parameter alpha is set to 0.51. Rerun
the script. What happens?
Reset alpha to 0.49, set plotit to 0 and maxlevel to 9. Rerun
the script. What can you say about the behaviour of the error norm?
Can you explain the behaviour, particularly in view of the nature of
the truncation error?
This scheme also has \( O(\Delta t, \Delta x^2) \) truncation error.
DOWNLOAD SOURCE
The driver script tdiff_1d_imp.m is configured so that as level is
varied, the ratio \(\Delta t/\Delta x\) (dtbydx) is held fixed.
EXERCISE
Run the script. What is the behaviour of the error norm in this case, and can you explain it?
Now study the function diff_id_imp and ensure that you understand how
it works, especially how the tridiagonal system is set up and solved.
Note that the linear system that needs to be solved is the
same at each time step, so we can define the sparse matrix \( A \)
outside of the time stepping loop. On the other hand, the right hand
side of the system changes from step to step (since it depends on
\( u^n \)), so its definition needs to be inside the loop.
EXERCISE
Working from diff_1d_imp.m and tdiff_1d_imp.m, implement a solution
of the diffusion equation using the Crank-Nicolson scheme:
\[ \frac{u^{n+1}_j - u^n_j}{\Delta t} = \frac 12 \left( \frac{u^{n+1}_{j+1}-2u^{n+1}_j+u^{n+1}_{j-1}}{\Delta x^2} + \frac{u^{n}_{j+1}-2u^{n}_j+u^{n}_{j-1}}{\Delta x^2} \right), \quad\quad j = 2,3,\ldots n_x -1 \] \[ u^{n+1}_1 = u^{n+1}_{n_x} = 0 \]
To get started, rewrite the above equation in the form \[ c_+ u^{n+1}_{j+1} + c_0 u^{n+1}_j + c_- u^{n+1}_{j-1} = f_j \] where \(c_+\), \(c_0\) and \(c_-\) are coefficients.
What can you say about the error norm as a function of discretization level for this scheme?