This post comes to serve as a complement to a previous one on the Runge-Kutta method from an algorithm: The 4th-order Runge-Kutta method for systems of ode in F90.
In that post, I presented a code in F90 to solve a system of 1st order ODE's. This is in agreement with was is usually presented in textbooks like: Numerical Methods for Engineers by Chapra, SC and Canale, RP and Differential Equations by Zill, DG. However, most of the differential equations found by scientists and engineers are not of 1st order!
This is quite enough when dealing with initial value problems but if boundary value problems appear some fix need to be introduced.
How do you handle a higher order ODE?
Easy: you need to split the equation into a set of 1st order ODE's.
This easlily done by means of one or several variable changes. From the point of view of the ODE's input to the program our problem is solved. However, another one arises: more conditions are needed and possibly unknown. This last subject is discussed in a section ahead.
Let us first tackle the first situation about the higher order of the ODE. Consider the following 5th-order ODE,
$2\dfrac{d^5y}{dx^5}+\dfrac{3}{5}\dfrac{d^3y}{dx^3}-\dfrac{dy}{dx}+9y=8x-1$ Eq. (01)
Equation (01) can be split by using the following changes of variables,
$y_1=y$
and for the first derivative of $y$
$\dfrac{dy_1}{dx}=\dfrac{dy}{dy}=y_2$and for the second derivative of $y$
$\dfrac{dy_2}{dx}=\dfrac{d^2y}{dx^2}=\dfrac{dy_1}{dx}=y_3$and for the third derivative of $y$
$\dfrac{dy_3}{dx}=\dfrac{d^3y}{dx^3}=\dfrac{dy_2}{dx}=y_4$and for the fourth derivative of $y$
$\dfrac{dy_4}{dx}=\dfrac{d^4y}{dx^4}=\dfrac{dy_3}{dx}=y_5$and for the fifth derivative of $y$
$\dfrac{dy_5}{dx}=\dfrac{d^5y}{dx^5}$
Next, Eq. (01) should now be written as the system of ODE,
$\dfrac{dy_1}{dx}=y_2$
$\dfrac{dy_2}{dx}=y_3$
$\dfrac{dy_3}{dx}=y_4$
$\dfrac{dy_4}{dx}=y_5$
$\dfrac{dy_5}{dx}=4x-\dfrac{1}{2} - \dfrac{3}{10}y_4 + \dfrac{1}{2}y_2 - \dfrac{9}{2}y_1$ Eqs. (02)
The system of ODE's in (02) are exactly the same posed in Eq. (01). Notice that in Eqs. (02) there are only first order ODE and which is more important: these are given explicitly, just as needed by the 4th-order Runge-Kutta method!
What about the conditions?
Let us consider a typical set of conditions for Eq. (01) such as,
$y=0$ at $x=0,1$
$\dfrac{dy}{dx}=0$ at $x=0,1$
$\dfrac{d^2y}{dx^2}=0$ at $x=1$
for example. If we translate these conditions into our new set of variables we get the following,
$y_1=0$ at $x=0,1$
$y_2=0$ at $x=0,1$
$y_4=0$ at $x=1$
Now, if we check against the set of ODE's in Eqs. (02) it is noticeable that conditions for $y_3$ and $y_5$ are required. What do we do then?
Well, this is one of the things our friends Chapra and Canale do not discuss, in to much detail, in their textbook. The fix to this difficulty is also known as the shooting method.
IVP vs BVP from the numerical perspective
Initial value and boundary value problems are quite different from the numerical perspective as you may have seen. BVP's are, usually, solved by the shooting method which is based on a linear relationship of the results for guesses of the unknown conditions so that a linear interpolation can be used (just as explained by Chapra and Canale).
However, the problem may become difficult to handle if complex numbers and functions appear. In some other cases, boundary conditions may be mixed (Robin type). Also, the ODE may be non linear so that different changes of variables may be needed as well.
In other words, a simple interpolation may not always work. Numerical solution of ODE is an art.
The eigenvalue problem
Here is another difficulty (if you thought that was all).
Sometimes, ODE are models that have to be adjusted for a given parameter so that you can come up with a solution of it. In this cases, you have this parameter (aka eignevalue) as part of the unknowns. What do you do in such a case?
It depends on the nature of every problem and of the obtain solutions. The solution must make physical sense, for example. In most cases, the conditions are set to fix value and the eigenvalue is used as the shooting parameter, which can be complex or real or integer.
Finally, when working with BVP's it is also common to use more robust algorithms and routines. For example, the shooting method is usually complemented by the Newton-Raphson method for complex variables.
- Up and running with Numerical Recipes in Fortran 90
- Bracketing roots (zbrac) - Numerical Recipes in F90
- An improvement for bracketing roots (zbrak) - Numerical recipes in F90
- The bisection method (rtbis) - Numerical Recipes in F90
- The false position method (rtfls) - Numerical recipes in F90
- The 4th-order Runge-Kutta method for systems of ode in F90
No comments:
Post a Comment