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 solve the system with the LU method and display the time that the LU algorithm take to solve the system using the timer Python module. Notice that, for \(N\) larger than about 340, the algorithm takes more than the maximum time allowed to run on the web (60 seconds) and therefore you won't get any output. On your own PC you can set \(N\) to higher values and, in general, it will be faster. For instance, on my desktop (with a i5-6500 CPU @ 3.20GHz) the time spent for \(N=300\) is about 8 seconds, while on the web it takes about 34 seconds. For \(N=400\) it takes about 20 seconds, while on the web it takes more than a minute, and for \(N=1000\) it takes about 5 minutes. It is quite interesting comparing these times with the ones achieved by the SciPy own solving algorithm in next section!