Section 4.7 Example: A Boundary Value Problem
As a concrete example, let us use the methods illustrated above to solve the following Ordinary Differential Equation:
\begin{equation*}
x''=x-4te^t,\;\;x(0)=0,\;x(1)=0,
\end{equation*}
whose exact solution is
\begin{equation*}
x(t)=t(1-t)e^t.
\end{equation*}
A way to solve numerically this ODE is to discretize the time interval \([0,1]\) into \(N\) segments of length \(h=1/N\) with endpoints \(t_0,t_1,\dots,t_N\text{.}\) We denote by \(x_0,x_1,\dots,x_N\) the values of the solution \(x(t)\) at the times \(t_i\text{.}\) Notice that \(x_0=0\) and \(x_N=0\) for each value of \(N\text{.}\) We discretize the second derivative analogously to how we did already for the first derivative in Section 2.2, namely we use the forward differencing approximation for the first derivative
\begin{equation*}
x'_i = x'(t_i) = \frac{x_{i+1}-x_i}{h}
\end{equation*}
and then the backward differencing for the second derivative:
\begin{equation*}
x''_i = x''(t_i) = \frac{x'_{i}-x'_{i-1}}{h} = \frac{x_{i+1}-2x_i+x_{i-1}}{h^2}.
\end{equation*}
Then the ODE above becomes a linear system in \(N-1\) unknowns:
\begin{equation*}
x_{i+1}-2x_i+x_{i-1} = h^2 (x_i - 4t_ie^{t_i}), \;\;i=1,\dots,N-1.
\end{equation*}
In matricial notation, the linear system looks as follows:
\begin{equation*}
\begin{pmatrix}
2+h^2& -1& 0 & & &\cr
-1&2+h^2 &-1&\ddots&\cr
0 &-1&\ddots&-1&0\cr
&\ddots&-1&2+h^2 &-1\cr
&&0&-1&2+h^2\cr
\end{pmatrix}
\begin{pmatrix}
x_1\cr x_2\cr\vdots\cr x_{N-1}\cr
\end{pmatrix}
=
4h^2
\begin{pmatrix}
t_1e^{t_1}\cr t_2e^{t_2}\cr\vdots\cr t_{N-1}e^{t_{N-1}}\cr
\end{pmatrix}
\end{equation*}
The code below solves numerically the system above with the three methods we covered in this chapter. The variable N is the number of subdivisions of the interval \([0,1]\) and the variable M is the maximum number of iterations we do with the iterative methods. Then we plot the exact solution (in blue), the numerical solution of the Jacobi method (in red), of the Gauss-Seidel method (in green) and of the SOR method (in cyan). Below this figure we also plot the error in the three methods with the relative colors. Notice that, for low values of N, the error in the iterative methods versus the number of iterations reaches soon a plateau, namely from that moment on the precision increases very slowly with the number of iterations. In order to increase the precision it is then needed to increase N. In the code below we add code to solve the system with the LU method. Since this method is not iterative, we plot the error as a horizontal line (in yellow). Note that in this case, when N grows "too much", the LU method gets heavy and slow and shortly after \(N=30\) the code won't be allowed to run in the web page because it takes too many resources (you can copy the code and run it in your MATLAB/Octave installation to try it with higher values of \(N\)).