Question 7.1: Use the method of lines to construct a numerical solution to...

The first step is the definition of the spatial discretization step
Δx =\frac{ x_r − x_l}{n − 1 }                                    (7.2.4)
At each node i, there is now a time-dependent variable T_i(t) that we seek to determine. Using a centered finite difference approximation of the spatial derivative for an interior node, the PDE in Eq. (7.2.1) becomes

\frac{dT_i}{dt}=\alpha \left[\frac{T_{i+1}(t) − 2T_i(t) + T_{i−1}(t)}{\Delta x^2} \right]                                            (7.2.5)

Note that this is now a system of coupled ODEs for the n variables T_i, since these discrete variables no longer depend on position. The initial condition at the interior nodes is the initial condition for the original PDE (7.2.3),

T_i(t = 0) = 2                                        (7.2.6)

The Dirichlet conditions imposed in Eqs. (7.2.2) suggest that the boundary nodes have a constant value for all times, i.e. T_1 = 0  and  T_n = 1. There are two ways that we can handle this situation. One option is to substitute these fixed values into the ODEs for the i = 2 and i = n − 1 nodes. As a result, Eq. (7.2.5) for node i = 2 becomes

\frac{dT_2}{dt}=\alpha \left[\frac{T_3(t)-2T_2(t)}{\Delta x^2} \right]                     (7.2.7)

and for node i = n − 1 we have

\frac{dT_{n-1}}{dt}=\alpha \left[\frac{1 − 2T_{n−1}(t) + T_{n−2}(t)}{\Delta x^2} \right]                           (7.2.8)

We then proceed by integrating the system of n − 2 ODEs embodied by Eqs. (7.2.5)– (7.2.8). Alternatively, and perhaps more conveniently, the following ODEs can be directly postulated for the boundary nodes,

\frac{dT_1}{dt} = 0,   T_1(0) = 0                                (7.2.9)

\frac{dT_n}{dt} = 0,  T_n(0) = 1                                 (7.2.10)

which together with the ODEs for the interior nodes lead to a system of n coupled ODEs in n unknown functions T_i(t). In compact form these will be

\frac{dT}{dt} = f(T),                      T(0) = T_0                            (7.2.11)

where T is the vector of the unknown variables

T=\begin{bmatrix} T_1 \\ T_2 \\ \vdots \\ T_n \end{bmatrix}                              (7.2.12)

f =\frac{ α}{Δx^2}  \begin{bmatrix} 0 \\ (T_3 − 2T_2 + T_1) \\ \vdots \\ (T_{i+1} − 2T_i + T_{i−1}) \\ \vdots \\ (T_n − 2T_{n−1} + T_{n−2})\\ 0 \end{bmatrix}                             (7.2.13)

and the initial condition vector is

T_0=\begin{bmatrix} 0 \\ 2 \\ \vdots \\ 2 \\ 1 \end{bmatrix}                              (7.2.14)

We can then integrate Eq. (7.2.11) using any of the methods we studied for systems of ordinary differential equations (e.g., RK4).