Question 7.3: Write a MATLAB program that uses implicit Euler and centered...
Write a MATLAB program that uses implicit Euler and centered finite differences to solve the diffusion–reaction equation
\frac{∂c}{∂t} = D \frac{∂^2c}{∂x^2} − Kc^2 (7.2.31)
subject to c(0, t) = 1, c(1, t) = 1, and c(x, 0) = 1.
Recall from our discussion in Section 4.6.2 that the steady state residual for a system of equations leads to a compact way to write the implicit Euler algorithm. It is thus useful to start with the steady state solution of this problem, i.e. the solution to the equation
0= D \frac{∂^2c}{∂x^2} − Kc^2 (7.2.32)
subject to the boundary conditions above. We first define n equidistant nodes in the interval [0, 1]. For the first node, from the boundary condition we have
c_1 = 1 (7.2.33)
For the interior nodes, using a centered finite difference approximation we get
D\left(\frac{c_{i+1} − 2c_i + c_{i−1}}{\Delta x^2} \right) − Kc^2_i = 0 (7.2.34)
For the last node, again from the boundary condition, we have
c_n = 1 (7.2.35)
These equations in residual form become
R^{ss}_1 = c_1 − 1 (7.2.36)
R^{ss}_i = \frac{D}{Δx^2} c_{i+1} − \frac{ 2D}{Δx^2} c_i − Kc^2_i +\frac{D}{Δx^2} c_{i−1} (7.2.37)
R^{ss}_n = c_n − 1 (7.2.38)
A MATLAB file to compute this residual is:
1 function R = getR(y,paramet)
2 dx = paramet(1);
3 D = paramet(2);
4 K = paramet(3);
5 n = length(y);
6 R = zeros(n,1);
7 R(1) = y(1) - 1;
8 for i = 2:n-1
9 R(i) = D/dx^2*(y(i+1) - 2*y(i) + y(i-1)) - K*y(i)^2;
10 end
11 R(n) = y(n) - 1;
As usual, we also need the steady state Jacobian. The non-zero Jacobian entries on the first row are
J^{ss}_{1,1} = 1 (7.2.39)
and on the last row are
J^{ss}_{n,n} = 1 (7.2.40)
For the interior rows, the non-zero entries are
J^{ss}_{i,i−1} = \frac{D}{Δx^2} (7.2.41)
J^{ss}_{i,i} = − \frac{2D}{Δx^2} − 2Kc_i (7.2.42)
J^{ss}_{i,i+1} = \frac{D}{Δx^2} (7.2.43)
The Jacobian is evaluated from the MATLAB function:
1 function J = getJ(y,paramet)
2 dx = paramet(1);
3 D = paramet(2);
4 K = paramet(3);
5 n = length(y);
6 J = zeros(n);
7 J(1,1) = 1;
8 for i = 2:n-1
9 J(i,i-1) = D/dx^2;
10 J(i,i) = -2*D/dx^2 - 2*K*y(i);
11 J(i,i+1) = D/dx^2;
12 end
13 J(n,n) = 1;
If we want to compute the steady state solution, we can use Newton–Raphson:
1 function [x,flag] = nonlinear_newton_raphson(x,tol,paramet)
2 % x = initial guess for x
3 % tol = criteria for convergence
4 % requires function R = getR(x)
5 % requires function J = getJ(x)
6 % flag = 0 if failed, 1 if converged
7 % paramet = [dx, D, K]
8 R = getR(x,paramet);
9 k = 0; %counter
10 kmax = 100; %max number of iterations allowed
11 while norm(R) > tol
12 J = getJ(x,paramet); %compute Jacobian
13 J = sparse(J);
14 del = -J\R; %Newton Raphson
15 x = x + del; %update x
16 k = k + 1; %update counter
17 R = getR(x,paramet);
18 if k > kmax
19 fprintf('Did not converge.\n')
20 break
21 end
22 end
23 if k > kmax || max(x) == Inf || min(x) == -Inf
24 flag = 0;
25 else
26 flag = 1;
27 end
Let’s move on to the unsteady problem. For each node, we derive the ODE that describes the evolution of the variable in time. For the first node this is
\frac{dc_1}{dt} = 0 (7.2.44)
For an interior node i, it is
\frac{dc_i}{dt} =\frac{D}{Δx^2} c_{i+1} − \frac{ 2D}{Δx^2} c_i − Kc^2_i +\frac{D}{Δx^2} c_{i−1} (7.2.45)
For the last node it is
\frac{dc_n}{dt} = 0 (7.2.46)
The initial condition is c_i(0) = 1 for all i.
The implicit Euler iterative scheme for propagating the solution in time takes the form
R_1 = c^{(k+1)}_1 − c^{(k)}_1 (7.2.47)
R_i= c^{(k+1)}_i − c^{(k)}_i - h\left[\frac{D}{Δx^2}(c^{(k+1)}_{i+1} − 2c^{(k+1)}_i + c^{(k+1)}_{i−1} ) − K(c^{(k+1)}_i )^2\right] (7.2.48)
R_n = c^{(k+1)}_n − c^{(k)}_n (7.2.49)
If you look at these equations, you will see that they look remarkably similar to the steady state solution. The residual for the interior nodes can also be written in the shorthand notation as
c^{(k+1)} − c^{(k)} − hR^{(k+1)}_{ss} = 0 (7.2.50)
where R^{(k+1)}_{ss} refers to the steady state residual. Likewise, we can reuse our Jacobian for the steady state problem since the partial derivatives of Eq. (7.2.50) with respect to c^{(k+1)}_i are simply
J = I − hJ^{(k+1)}_{ss} (7.2.51)
The easiest way to combine the time integration with what we wrote for the steady state is to modify our program for implicit Euler for a system of equations (Program 4.5) to mesh with the functions we already wrote for getR and getJ:
1 function [xout,yout] = ivp_implicit_euler_sys(x0,y0,h,n_out,i_out,p)
2 % x0 = initial value for x
3 % y0 = initial value for y
4 % h = step size
5 % i_out = frequency to output x and y to [xout, yout]
6 % n_out = number of output points
7 % requires the function getR(y,yold,h) and getJ(y,h)
8 % p = [dx, D, K]
9 xout = zeros(n_out+1,1); xout(1) = x0;
10 yout = zeros(n_out+1,length(y0)); yout(1,:) = y0;
11 x = x0; y = y0; tol = 1e-8;
12 for j = 2:n_out+1
13 for i = 1:i_out
14 yold = y;
15 R=y - yold - h*getR(y,p);
16 while norm(R) > tol
17 J = eye(length(y)) - h*getJ(y,p);
18 del = -J\R;
19 y = y + del;
20 R=y - yold - h*getR(y,p);
21 end
22 x = x + h;
23 end
24 xout(j) = x;
25 yout(j,:) = y;
26 end
Notice that Eq. (7.2.50) is implemented on lines 15 and 20, while Eq. (7.2.51) is implemented on line 17. Figure 7.3 shows the integration of this problem for a Thiele modulus of \sqrt{2} .
