Question 8.5: FINITE-DIFFERENCE METHOD, LINEAR BVP WITH MIXED BOUNDARY CON...
FINITE-DIFFERENCE METHOD, LINEAR
BVP WITH MIXED BOUNDARY CONDITIONS
Consider the BVP
t \ddot{w}+\dot{w}+2 t=0, \quad w(1)=1, \quad \dot{w}(2.5)=-1
Solve by the finite-difference method and \Delta t=0.25. Use central-difference approximations for all derivatives. Compare the results with the exact solution values at the grid points. Also confirm that the boundary condition at the right end is satisfied by the computed solution.
The interval is 1.5 in length and \Delta t=0.25, hence N=1.5 / 0.25=6, and there are five interior grid points. Replacing the first and second derivatives with central-difference formulas yields
t_{i} \frac{w_{i-1}-2 w_{i}+w_{i+1}}{\Delta t^{2}}+\frac{w_{i+1}-w_{i-1}}{2 \Delta t}+2 t_{i}=0
Simplify the above and apply at interior grid points so that
\left(2 t_{i}-\Delta t\right) w_{i-1}-4 t_{i} w_{i}+\left(2 t_{i}+\Delta t\right) w_{i+1}=-4 t_{i} \Delta t^{2}, \quad i=2,3,4,5,6
Consequently,
\begin{array}{ll} i=2 & \left(2 t_{2}-\Delta t\right) w_{1}-4 t_{2} w_{2}+\left(2 t_{2}+\Delta t\right) w_{3}=-4 t_{2} \Delta t^{2} \\ i=3 & \left(2 t_{3}-\Delta t\right) w_{2}-4 t_{3} w_{3}+\left(2 t_{3}+\Delta t\right) w_{4}=-4 t_{3} \Delta t^{2} \\ i=4 & \left(2 t_{4}-\Delta t\right) w_{3}-4 t_{4} w_{4}+\left(2 t_{4}+\Delta t\right) w_{5}=-4 t_{4} \Delta t^{2} \\ i=5 & \left(2 t_{5}-\Delta t\right) w_{4}-4 t_{5} w_{5}+\left(2 t_{5}+\Delta t\right) w_{6}=-4 t_{5} \Delta t^{2} \\ i=6 & \left(2 t_{6}-\Delta t\right) w_{5}-4 t_{6} w_{6}+\left(2 t_{6}+\Delta t\right) w_{7}=-4 t_{6} \Delta t^{2} \end{array} (8.7)
w_{1}=w(1)=1 is provided by the boundary condition at the left endpoint. At the right endpoint, however, w_{7}=w(2.5) is not directly available, but \dot{w}(2.5) is. To approximate \dot{w} at the right end, we will use a one-sided, three-point backward difference formula so that the values at the previous points are utilized. Note that this has second-order accuracy (see Table 6.3, Section 6.2), which is in line with the central-difference formulas used in the earlier stages.
\left.\dot{w}\right|_{t_{i}} \cong \frac{w_{i-2}-4 w_{i-1}+3 w_{i}}{2 \Delta t} (8.8)
Applying Equation 8.8 at the right end (i=7),
\frac{w_{5}-4 w_{6}+3 w_{7}}{2 \Delta t}=-1 \stackrel{\text { Solve for } w_{7}}{\Rightarrow} w_{7}=\frac{1}{3}\left(-2 \Delta t+4 w_{6}-w_{5}\right) (8.9)
Substitution into Equation 8.7 for w_{7}, and expressing the system in matrix form, yields
\left[\begin{array}{ccccc} -4 t_{2} & 2 t_{2}+\Delta t & 0 & 0 & 0 \\ 2 t_{3}-\Delta t & -4 t_{3} & 2 t_{3}+\Delta t & 0 & 0 \\ 0 & 2 t_{4}-\Delta t & -4 t_{4} & 2 t_{4}+\Delta t & 0 \\ 0 & 0 & 2 t_{5}-\Delta t & -4 t_{5} & 2 t_{5}+\Delta t \\ 0 & 0 & 0 & \frac{4}{3}\left(t_{6}-\Delta t\right) & -\frac{4}{3}\left(t_{6}-\Delta t\right) \end{array}\right]\left\{\begin{array}{c} w_{2} \\ w_{3} \\ w_{4} \\ w_{5} \\ w_{6} \end{array}\right\}=\left\{\begin{array}{c} -4 t_{2} \Delta t^{2}-\left(2 t_{2}-\Delta t\right) w_{1} \\ -4 t_{3} \Delta t^{2} \\ -4 t_{4} \Delta t^{2} \\ -4 t_{5} \Delta t^{2} \\ -4 t_{6} \Delta t^{2}+\frac{2}{3}\left(2 t_{6}+\Delta t\right) \Delta t \end{array}\right\} (8.10)
The coefficient matrix is once again tridiagonal. This system is solved for the solution vector to obtain the approximate values at the interior grid points. The process is completed by attaching w_{1}, which is available, and w_{7}, which is found by Equation 8.9 , to the vector of the interior grid values. The following MATLAB script performs these tasks:
TABLE 8.3
Comparison of Exact and Computed Values at Grid Points (Example 8.5)
| Solution | t = 1 | t = 1.25 | t = 1.5 | t = 1.75 | t = 2 | t = 2.25 | t = 2.5 |
| Computed | 1.0000 | 1.5599 | 1.9044 | 2.0804 | 2.1164 | 2.0304 | 1.8351 |
| Exact | 1.0000 | 1.5555 | 1.8955 | 2.0673 | 2.0993 | 2.0097 | 1.8111 |
a = 1; b = 2.5; dt = 0.25;
N = (b-a)/dt; t = a:dt:b;
w1 = 1; % BC at the left end
A = zeros(N-1,N-1); b = zeros(N-1,1); % Pre-allocate
b(1) = -4*t(2)*dt^2-(2*t(2)-dt)*w1;
b(N-1) = -4*t(N)*dt^2+(2/3)*(2*t(N)+dt)*dt;
for i = 1:N-2,
A(i,i+1) = 2*t(i+1)+dt;
end
for i = 1:N-3,
A(i+1,i) = 2*t(i+2)-dt;
b(i+1) = -4*t(i+2)*dt^2;
end
A(N-1,N-2) = (4/3)*(t(N)-dt);
for i = 1:N-2,
A(i,i) = -4*t(i+1);
end
A(N-1,N-1) = -(4/3)*(t(N)-dt);
w_interior = ThomasMethod(A,b);
w = [w1 w_interior'];
w(N+1) = (1/3)*(-2*dt+4*w(N)-w(N-1));
>> w
w =
1.0000 1.5599 1.9044 2.0804 2.1164 2.0304 1.8351
% Exact solution
>> w_ex = matlabFunction(dsolve('t*D2w+Dw+2*t=0','w(1)=1','Dw(2.5)=-1'));
>> w_e = w_ex(t)
w_e =
1.0000 1.5555 1.8955 2.0673 2.0993 2.0097 1.8111
Table 8.3 summarizes the computed and exact values at the grid points.
In order to confirm that the solution obtained here satisfies the boundary condition at the right end, we run the same code with t=0.0125 to generate a smooth solution curve. The result is shown in Figure 8.4.
