Question 8.4: FINITE-DIFFERENCE METHOD, NONLINEAR BVP Consider the BVP in ...
FINITE-DIFFERENCE METHOD, NONLINEAR BVP
Consider the BVP in Example 8.2:
T_{x x}-\beta T^{4}-\alpha T+\gamma=0, \quad T(0)=500, \quad T(0.2)=350
where \alpha=20, \beta=10^{-8}, and \gamma=5 \times 10^{3} with all parameters in consistent physical units. Solve by the finite-difference method using central-difference formulas and h=\Delta x=0.04.
The interval is 0.2 in length and \Delta x=0.04, hence N=0.2 / 0.04=5, and there are N-1=4 interior grid points. At each interior grid point, the second derivative is replaced with a central-difference formula to obtain
\frac{T_{i-1}-2 T_{i}+T_{i+1}}{\Delta x^{2}}-\beta T_{i}^{4}-\alpha T_{i}+\gamma=0
Simplifying the above, we have
T_{i-1}-\left(2+\alpha \Delta x^{2}\right) T_{i}+T_{i+1}-\beta \Delta x^{2} T_{i}^{4}+\gamma \Delta x^{2}=0, \quad i=2,3,4,5 (8.5)
Applying Equation 8.5, keeping in mind that T_{1}=500 and T_{6}=350, we find
\begin{gathered} 508-2.0320 T_{2}-1.6 \times 10^{-11} T_{2}^{4}+T_{3}=0 \\ T_{2}-2.0320 T_{3}-1.6 \times 10^{-11} T_{3}^{4}+T_{4}+8=0 \\ T_{3}-2.0320 T_{4}-1.6 \times 10^{-11} T_{4}^{4}+T_{5}+8=0 \\ 358+T_{4}-2.0320 T_{5}-1.6 \times 10^{-11} T_{5}^{4}=0 \end{gathered}
We will solve this system of nonlinear equations using Newton’s method for systems; see Chapter 4 . To that end, we first conveniently express the above system as
\begin{aligned} & f_{1}\left(T_{2}, T_{3}, T_{4}, T_{5}\right)=0 \\ & f_{2}\left(T_{2}, T_{3}, T_{4}, T_{5}\right)=0 \\ & f_{3}\left(T_{2}, T_{3}, T_{4}, T_{5}\right)=0 \\ & f_{4}\left(T_{2}, T_{3}, T_{4}, T_{5}\right)=0 \end{aligned}
Newton’s method requires that in each step we solve a linear system in the form
\left[\begin{array}{c}\frac{\partial f_1}{\partial T_2} & \cdot & \cdot & \frac{\partial f_1}{\partial T_5} \\ \cdot & & & \cdot \\ \cdot & & & \cdot \\ \frac{\partial f_4}{\partial T_2} & \cdot & \cdot & \frac{\partial f_4}{\partial T_5} \end{array}\right] \left\{\begin{array}{c} \Delta T_2 \\ \cdot \\ \cdot \\ \Delta T_5 \end{array}\right\}=\left\{\begin{array}{c} -f_1 \\ \cdot \\ \cdot \\ -f_4 \end{array}\right\} (8.6)
It is readily observed that the coefficient matrix in Equation 8.6 is tridiagonal, thus the system in Equation 8.6 is most efficiently solved using the Thomas method.
Proceeding as in Example 4.20, the following script implements Newton’s method to find the solution estimates at the four interior grid points. The initial values are arbitrarily selected as 400 . The maximum of number of iterations, as well as the tolerance are chosen as in Example 4.20.
TABLE 8.2
Comparison of Results: Shooting Method, Finite Difference (Example 8.4)
| x | Finite Difference | Shooting Method |
| 0 | 500 | 500 |
| 0.04 | 457.6786 | 457.6375 |
| 0.08 | 422.7049 | 422.6493 |
| 0.12 | 393.7686 | 393.7178 |
| 0.16 | 369.8176 | 369.7862 |
| 0.2 | 350 | 350 |
syms T2 T3 T4 T5
f1 = 508-2.0320*T2-1.6e-11*T2^4+T3; f2 = T2-2.0320*T3-1.6e-11*T3^4+T4+8;
f3 = T3-2.0320*T4-1.6e-11*T4^4+T5+8; f4 = 358+T4-2.0320*T5-1.6e-11*T5^4;
f = [f1;f2;f3;f4];
J = matlabFunction(jacobian(f,[T2,T3,T4,T5]));
F = matlabFunction(f);
tol = 1e-4; kmax = 20; T(1,:) = 400*ones(1,4);
for k = 1:kmax,
A = J(T(k,1),T(k,2),T(k,3),T(k,4));
b = -F(T(k,1),T(k,2),T(k,3),T(k,4));
if abs(b(1)) < tol && abs(b(2)) < tol && abs(b(3)) < tol && abs(b(4))
< tol,
root = T(k,:);
return
end
if det(A) == 0,
break
end
delT = ThomasMethod(A,b); delT = delT';
T(k+1,:) = T(k,:) + delT;
if norm(delT) < tol,
root = T(k+1,:);
break
end
end
Execution of this script shows that convergence is achieved after two iterations:
>> T
T =
400.0000 400.0000 400.0000 400.0000
457.7283 422.7499 393.8024 369.8407
457.6786 422.7049 393.7686 369.8176
Comparison with the Shooting Method
In order to compare these solution estimates with those obtained by the shooting method (Example 8.2), we adjust the step size used in RK4 to 0.04. A summary of the results is shown in Table 8.2, where it is observed that the estimates provided by the two techniques closely follow one another.