1-D Linear ODEs of 1st and 2nd order
The general forms of 1st and 2nd order ODEs are as follows:
\begin{equation} \begin{aligned} & U' + aU = f, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \in [\alpha,\beta], \ \ \text{plus 1 boundary condition} \\ & U'' + aU' + bU = f, \ \ \ x \in [\alpha,\beta], \ \ \text{plus 2 boundary conditions} \end{aligned} \end{equation}
in which $ a $ and $ b $ can be functions of $ x $. To solve these types of differential equations numerically, the simplest way is to use the finite difference approximation. Using Taylor series, different schemes could be easily derived. Below the Taylor expansion series are presented for three adjacent nodes:
\begin{equation} \begin{aligned} & U_{j-1} = U_j + U'_j (x_{j-1}-x_j) + U''_j \frac{(x_{j-1}-x_j)^2}{2!} + U^{(3)}_j \frac{(x_{j-1}-x_j)^3}{3!} + \cdots \\ \\ & U_j = U_j \\ & U_{j+1} = U_j + U'_j (x_{j+1}-x_j) + U''_j \frac{(x_{j+1}-x_j)^2}{2!} + U^{(3)}_j \frac{(x_{j+1}-x_j)^3}{3!} + \cdots \end{aligned} \end{equation}
In order to simplify the formulas, it holds:
\begin{equation} \Delta x_n = x_n - x_{n-1} \end{equation}
1st order derivatives
For the 1st order ODEs, three schemes will be employed; backward, centered and forward schemes. To derive the backward scheme, the first and second expressions mentioned above are linearly combined:
\begin{equation} A \ U_{n-1} + B \ U_n = A \left( U_n - U'_n \Delta x_n + U''_n \frac{(\Delta x_n)^2}{2!} + \cdots \right) + B \ U_n, \end{equation}
where $ A $ and $ B $ are suitable constants.
To assure the consistency of the scheme, the term of order $ O\left((\Delta x_n)^0\right) $ on the right hand side must be erased ($ A + B = 0 $); furthermore, to directly express the first order derivative, we can set $ B = 1/{\Delta x_n} $. Then we have:
\begin{equation} U' \vert _{x = x_n} = \frac{U_n - U_{n-1}}{\Delta x_n} + O(\Delta x_n). \end{equation}
Proceeding in the same way, the forward scheme can be obtained using the second and third relations:
\begin{equation} U' \vert _{x = x_n} = \frac{U_{n+1} - U_n}{\Delta x_{n+1}} + O(\Delta x_{n+1}); \end{equation}
And for the centered scheme, one should combine first and third relations:
\begin{equation} U' \vert _{x = x_n} = \frac{U_{n+1} - U_{n-1}}{\Delta x_{n+1} + \Delta x_n} + O(\Delta x_{n+1} + \Delta x_n); \end{equation}
In this case if an evenly spaced grid is used, the scheme order increases to two:
\begin{equation} U' \vert _{x = x_n} = \frac{U_{n+1} - U_{n-1}}{2 \Delta x} + O\left((\Delta x)^2\right). \end{equation}
2nd order derivatives
For the 2nd order ODEs, for efficiency reasons in terms of computational costs and convergence speed, just the centered scheme will be employed. To derive this scheme, all three expressions (Taylor expansions) are linearly combined:
\begin{equation} A \ U_{n-1} + B \ U_n + C \ U_{n+1} = A \left( U_n - U'_n \Delta x_n + U''_n \frac{(\Delta x_n)^2}{2!} + \cdots \right) + \\ + B \ U_n + C \left( U_n + U'_n \Delta x_{n+1} + U''_n \frac{(\Delta x_{n+1})^2}{2!} + \cdots \right), \end{equation}
where $ A $, $ B $ and $ C $ are suitable constants.
To assure the consistency of the scheme, the terms of orders $ O\left((\Delta x)^0\right) $ and $ O\left((\Delta x)^1\right) $ on the right hand side must be erased ($ A + B + C = 0 $, $ -A\Delta x_n + C\Delta x_{n+1} = 0 $). Furthermore, to directly express the second order derivative, we can set $ A \Delta x_n^2 + C \Delta x_{n+1}^2 = 2 $. Then we have:
\begin{equation} U'' \vert _{x = x_n} = \frac{\frac{U_{n+1} - U_n}{\Delta x_{n+1}} - \frac{U_n - U_{n-1}}{\Delta x_n}}{\frac{\Delta x_{n+1} + \Delta x_n}{2}} + O(\Delta x_{n+1} + \Delta x_n); \end{equation}
In this case if an evenly spaced grid is used, the scheme order increases to two:
\begin{equation} U'' \vert _{x = x_n} = \frac{U_{n+1} - 2 U_n + U_{n-1}}{\Delta x^2} + O\left((\Delta x)^2\right). \end{equation}
Discretizing the 1st order ODEs
Considering the general form of the 1st order linear ODEs:
\begin{equation} U' + aU = f, \ \ \ \ \ \ \ \ \ \ \ \ \ \ x \in [\alpha,\beta], \ \ \text{plus 1 boundary condition} \\ \end{equation}
its discretization by means of finite difference backward scheme is:
\begin{equation} \frac{U_n - U_{n-1}}{\Delta x_n} + aU_n + O(\Delta x_n) = f_n. \end{equation}
Here for simplicity, $ a $ is assumed to be constant. Also, for forward and centred schemes we have:
\begin{equation} \frac{U_{n+1} - U_n}{\Delta x_{n+1}} + aU_n + O(\Delta x_{n+1}) = f_n, \end{equation}
\begin{equation} \frac{U_{n+1} - U_{n-1}}{\Delta x_{n+1} + \Delta x_n} + aU_n + O(\Delta x_{n+1} + \Delta x_n) = f_n. \end{equation}
In order to check the implementation of the scheme, it is possible to consider a suitable first order ODE with a known analytical solution:
\begin{equation} U' + aU = \cos x + a \sin x, \ \ \ [0, \pi]. \end{equation}
When a homogeneous boundary condition is applied on the left, the exact solution is:
\begin{equation} U = \sin x, \end{equation}
and in this case, the convergence rate of the backward scheme could be verified.
Similarly, when the homogeneous boundary condition is applied on the right, we can test the convergence rate of the forward scheme.
Finally, for the centered scheme we can apply the boundary condition either on left or right, locally changing the scheme to backward or forward respectively. The convergence rate of the centered scheme is presented as follows.
The relevant Matlab codes can be downloaded using the following link:
Matlab Code: 1st_Order_Eq.tar.gz
Discretizing the 2nd order ODEs
Considering the general form of the 2nd order linear ODEs:
\begin{equation} U'' + aU' + bU = f, \ \ \ x \in [\alpha,\beta], \ \ \text{plus 2 boundary conditions} \\ \end{equation}
its discretization by means of finite difference centred scheme is:
\begin{equation} \frac{\frac{U_{n+1} - U_n}{\Delta x_{n+1}} - \frac{U_n - U_{n-1}}{\Delta x_n}}{\frac{\Delta x_{n+1} + \Delta x_n}{2}} + a \frac{U_{n+1} - U_{n-1}}{\Delta x_{n+1} + \Delta x_n} + bU_n + O(\Delta x_{n+1} + \Delta x_n) = f_n. \end{equation}
For the first order derivative, if instead of using a centred scheme it is implemented a forward or a backward scheme, we have:
\begin{equation} \frac{\frac{U_{n+1} - U_n}{\Delta x_{n+1}} - \frac{U_n - U_{n-1}}{\Delta x_n}}{\frac{\Delta x_{n+1} + \Delta x_n}{2}} + a \frac{U_{n+1} - U_n}{\Delta x_{n+1}} + bU_n + O(\Delta x_{n+1} + \Delta x_n) = f_n. \end{equation}
\begin{equation} \frac{\frac{U_{n+1} - U_n}{\Delta x_{n+1}} - \frac{U_n - U_{n-1}}{\Delta x_n}}{\frac{\Delta x_{n+1} + \Delta x_n}{2}} + a \frac{U_n - U_{n-1}}{\Delta x_n} + bU_n + O(\Delta x_{n+1} + \Delta x_n) = f_n. \end{equation}
As said in the previous sections, just using centred schemes for both derivatives is possible to increase the rate of convergence to order two evenly distributing the grid nodes:
\begin{equation} \frac{U_{n+1} - 2 U_n + U_{n-1}}{\Delta x^2} + a \frac{U_{n+1} - U_{n-1}}{2 \Delta x} + bU_n + O\left((\Delta x)^2\right) = f_n. \end{equation}
In order to check the implementation of the scheme, it is possible to consider a suitable second order ODE with a known analytical solution:
\begin{equation} U'' + aU' + bU = -\sin x + a \cos x + b \sin x, \ \ \ [0, 3 \pi/2]. \end{equation}
When Dirichlet's boundary conditions are suitably applied ($U(0) = 0 \ \ , \ \ U(3 \pi/2) = -1$), the exact solution is:
\begin{equation} U = \sin x, \end{equation}
and in this case, the convergence rate of the centred scheme on an evenly spaced grid could be verified.
When Neumann's boundary condition is applied on the left or on the right, discretizing is required. For simplicity it is possible to choose a FwdFD scheme for the left edge and a BwdFD scheme for the right one, but, whenever a first order scheme is chosen, it can affect the rate of convergence. Two examples are presented below:
The relevant Matlab codes can be downloaded using the following link:
Matlab Codes Bank
It is great for the beginner, thanks for sharing.
ReplyDelete