Question 10.6: FINITE-DIFFERENCE METHOD; 1D HEAT EQUATION Consider a latera...
FINITE-DIFFERENCE METHOD; 1D HEAT EQUATION
Consider a laterally insulated wire of length L=1 and \alpha=0.5, whose ends are kept at zero temperature, subjected to initial temperature f(x)=10 \sin \pi x. Compute the approximate values of temperature u(x, t), 0 \leq x \leq 1,0 \leq t \leq 0.5, at mesh points generated by h=0.25 and k=0.1. All parameter values are in consistent physical units. The exact solution is given as
u(x, t)=(10 \sin \pi x) e^{-0.25 \pi^{2} t}
We first calculate
r=\frac{k \alpha^{2}}{h^{2}}=\frac{(0.1)(0.5)^{2}}{(0.25)^{2}}=0.4<\frac{1}{2}
which signifies stability and convergence for the finite-difference method described by Equation 10.28.
u_{i, j+1}=(1-2 r) u_{i j}+r\left(u_{i-1, j}+u_{i+1, j}\right), \quad r=\frac{k \alpha^2}{h^2} (10.28)
With r=0.4, Equation 10.28 reduces to
u_{i, j+1}=0.2 u_{i j}+0.4\left(u_{i-1, j}+u_{i+1, j}\right) (10.30)
This will generate approximate solutions at the mesh points marked in Figure 10.11. The first application of Equation 10.30 is at the (1,0) position so that
u_{11}=0.2 u_{10}+0.4\left(u_{00}+u_{20}\right) \underset{\substack{\text { boundary } \\ \text { values }}}{\stackrel{\text { Using }}{=}} 0.2(7.0711)+0.4(0+10)=5.4142
It is next applied at (2,0) to find u_{21}, and at (3,0) to find u_{31}. This way, the values at the three mesh points along the time row j=1 become available. Subsequently,
Equation 10.30 is applied at the three mesh points on the current row (j=1) to find the values at the next time level, j=2, and so on. Continuing with this strategy, we will generate the approximate values shown in Figure 10.11. It is interesting to note that at the first time row, we are experiencing relative errors of around 2 \%, but this grows to about 9.63 \% by the time we arrive at the fifth time level.
The numerical results obtained in this manner can be readily verified by executing the user-defined function Heat1DFD:
>> t = 0:0.1:0.5; x = 0:0.25:1;
>> u0 = 10.*sin(pi*x);
>> u = Heat1DFD(t,x,u0,0.5)
u =
0 1.8610 2.6318 1.8610 0
0 2.4304 3.4372 2.4304 0
0 3.1742 4.4890 3.1742 0
0 4.1456 5.8627 4.1456 0
0 5.4142 7.6569 5.4142 0
0 7.0711 10.0000 7.0711 0
As stated earlier, the output resembles the gridded region.
