Here is, perhaps, the simplest case about a boundary value problem solved by a shooting method. You are warned that there are several shooting methods and there are even chapters and books devoted to these.
According to the nature of the system of ordinary differential equations and its corresponding boundary conditions some sort of shooting method can be applied/adpated.
Contents
1 - ODEs and its conditions
2 - Description of the numerical solution
3 - Algorithm to work with the shooting method
4 - The code
4.1 - Main program
4.2 - Subroutine to load initial vaues
4.3 - Subroutine for 4th-order Runge-Kutta method equations
4.4 - Subroutine for the evaluation of derivatives
5 - Results and comparison
6 - Possible improvements of the program
2 - Description of the numerical solution
3 - Algorithm to work with the shooting method
4 - The code
4.1 - Main program
4.2 - Subroutine to load initial vaues
4.3 - Subroutine for 4th-order Runge-Kutta method equations
4.4 - Subroutine for the evaluation of derivatives
5 - Results and comparison
6 - Possible improvements of the program
________________
1 - ODES and its conditions
Go to the topSince the idea is to implement a F90 program to solve a boundary value problem I have search for an example so that checks on the numerical solutions can be done. Fortunately, the book of Numerical Methods for Engineers by SC Chapra and RP Canale gives such an example in Chapter 27. However, that book does not give the algorithm nor the step-by-step process but the numerical results.
Therefore, the problem consists of the single 2nd-order differential equation,
subject to the following boundary conditions
at $x=0$, $T=40$
at $x=10$, $T=200$ Eq. (02)
where $T_A=20$ and $h=0.01$. This equation and its boundary conditions have an interesting physical background on the heating of a metal bar. Here, I shall focus on the mathematical aspects of it.
Since the variables in use are $y$ and $x$ in previous posts and supplied codes I shall switch to these so that the ODEs becomes,
$\dfrac{d^2 y(x)}{dx^2}=-0.01\left[ 20 - y(x) \right]$ Eq. (03)
subject to,
at $x=0$, $y=40$
at $x=10$, $y=200$ Eq. (04)
ODE Eq. (03) can be solved with the 4th-order Runge-Kutta (RK4) method and the boundary conditions can be satisfied with a shooting method. This is done as follows.
Equation (03) must be reduced to a set of 1st-order ODEs so that these can be suppplied to the RK4 program. However, since these program only accepts initial values guesses at $x=0$ should introduced (aka shoots) looking for the boundary condition at $x=10$ to be satisfied.
At this point one should wonder what a suitable guess could be. For short, there is no clue and we should proceed by trial and error.
If the change of variable,
$y=y_1$ Eq. (05)
is introduced, then the derivatives in Eq. (03) can be determined to be,
$\dfrac{dy}{dx}=\dfrac{dy_1}{dx}=y_2$
$\dfrac{d^2y}{dx^2}=\dfrac{d^2y_1}{dx^2}=\dfrac{dy_2}{dx}$ Eq. (06)
Next, Eq. (03) can be written as the following set of 1st-order ODE,
$\dfrac{dy_1}{dx}=y_2$
$\dfrac{dy_2}{dx}=-0.01\left[ 20 - y(x) \right]$ Eq. (07)
$\dfrac{dy_2}{dx}=-0.01\left[ 20 - y(x) \right]$ Eq. (07)
while the boundary conditions are translated as,
at $x=0$, $y_1=40$
at $x=10$, $y_1=200$ Eq. (08)
Equations (07-08) are now the problem at hand. A quick look to Eqns. (07-08) reveals that, at least, an initial condition for $y_2$ is missing. This why the shooting method is required since it must be determined along with the ODEs equations.
According to different book authors the shooting method is about supplying a suitable guess for the unknown initial value (or values since there may be more than one unknown!) so that the condition(s) at the other boundary are satisfied.
Of course, giving a guess value is easy but if you are not close to the exact value trial and error may consume too much time so that a more systematic way of doing is required. For this very simple case: an ordinary linear differential equation subject the shooting method is very simple:
Since the original ODE Eq. (03) is linear, the values of the $y(10)$ at two shootings are linearly related too. Therefore, a linear interpolation can be used to determine the exact value of the unknown function at $x=0$.
As you may imagine, the structure of the program should be changed since for this case it is the shooting method what is at the heart of the program because nothing can be calculated unless a guess for the unknown condition is supplied. Below is presented a flow diagram as a sort of algorithm,
 |
| Fig. 1 Flow diagram for the case a shooting method based on a trial and error approach to solve a boundary value problem |
For the present problem, the way in which the guess for the unknown condition from the two shootings i determined relies on the simple formula,
$y2_i(0)=10+\dfrac{y2_{i+1}(xi) - y2_{i-1}(xi)}{y1_{i+1}(xf) - y1_{i-1}(xf)}\left[ y1_{i}(xf)-y1_{i-1}(xf) \right]$ Eq. (09)
which can be easily found if the slope of the line connecting the 2 points, corresponding to the shootings, is considered and forcing a third point including the desired boundary condition to be satisfied be part of that line.
4 - The code
Go to the topThe program was structured following the idea of the flow diagram above. It cosists of a main program and a set of subroutines:
- load for loading the known "initial" conditions
- integrator for advancing steps
- RK4 to apply the 4th-order Runge-Kutta method formulas
- derivs for evaluation of the set of 1st-order ODE
Below is given the code so that it can readily implemented. I also use a module file to introduce the single and double precision kind parameters as well as f or integers. You may assemble this last module with,
MODULE nrtypes
!===============================================================================
!Symbolic names for kind types of 4-, 2-, and 1-byte integers:
!===============================================================================
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
!
!===============================================================================
!Symbolic names for kind types of single- and double-precision reals:
!===============================================================================
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
END MODULE nrtypes
The subroutines were also placed in another module file too. In this case, the module is called: MySubroutines.f90.
4.1 - Main program
Go to the topThe main program contains the initialization of some key parameters, callings to subroutines and printing of results. Her it is,
program ODEstestShooting
!
! This program is intended for testing the solving of a set of ODEs
! in an already solved example in the book of Chapra and Canale by ! means of a 4th order Runge-Kutta method.
!
! It only deals with the equations as an initial value problem with ! a single shooting.
!
use nrtypes
use MySubroutines
!
implicit none
integer(I4B), parameter:: n = 2 !Number of ode
integer(I4B), parameter:: m = 5 !Number of intervals from xi to xf
integer(I4B), parameter:: nn = 1 !Number of known boundary conditions
integer(I4B):: mc !Dummy counter for intervals from 0 to m
integer(I4B):: i, sc !Counter
real(SP):: x !Independent (x)
real(SP), dimension(m):: y !and independent (y) variables
real(SP), dimension(n):: yi !Initial values
real(SP), dimension(2):: y2g !Guessing values for unknwon y2 at xi
real(SP), dimension(2):: y1g !Obtain value from RK4 from the !shootings
real(SP):: xi, xf !Final values
real(SP):: dx !Step
real(SP):: xend
real(SP):: h !Dummy step
real(SP):: y1f ! Value for y(1) at xf
!
!Initialization of the variables and parameters
!Data of ode system
xi = 0.0_SP ! This is the initial value
xf = 10.0_SP ! This is the end of the integration range
y1f = 200.0_SP ! This is the value of y1(xf). The shootings aim to !this value
!Step size calculation
dx = (xf-xi)/m
!
x = xi
!
! The known boundary conditions at xi (initial conditions in fact) ! are supplied by a subroutine
call Load(yi)
!
do i = 1,nn
y(i) = yi(i)
!
end do
!
! Print a message for the user to supply two bracketing guesses for ! the unkown
! boundary condition: y2 in this case
PRINT*, 'Please, supply two guesses for shootings using y2'
read(*,*) y2g(1), y2g(2)
!
! Printing results
PRINT*, ''
PRINT*, '4th order Runge-Kutta pogram'
PRINT*, 'Uses shooting method'
!
PRINT*, '============RESULTS============'
PRINT*, ' x T'
PRINT*, '==============================='
PRINT*, x, y(1)
!
! Loop for shootings
do sc =1,2
call Load(yi)
do i = 1,nn
y(i) = yi(i)
end do
!
y(2) = y2g(sc)
x=xi
do mc = 1,m
xend = x + dx
if (xend > xf) then
xend = xf
endif
h = dx
call integrator(x,y,n,h,xend)
y1g(sc) = y(1)
end do
end do
!
! Next, check to see if the knwon condition y1(xf) is between the ! results
! from the two shootings
if (y1g(1)<y1f .and. y1g(2)>y1f) then
y(2) = y2g(1) + ( y2g(2) - y2g(1) )*( y1f - y1g(1) ) &
/( y1g(2) - y1g(1) )
else
STOP 'y1(xf) in not in the range. Give another gueses'
end if
!
! If y1(xf) has been correctly bracketed then the new y(2) can be ! used to
! calculate the correct numerical solution to the ODE
!
call Load(yi)
do i = 1,nn
y(i) = yi(i)
end do
x=xi
do mc = 1,m
xend = x + dx
if (xend > xf) then
xend = xf
endif
h = dx
call integrator(x,y,n,h,xend)
PRINT*, x, y(1)
if (x >= xf) exit
end do
!
PRINT*, '==============================='
!
end program
As you can see in the code, a conditional was introduced to check for correct bracketing of the aim boundary condition $y_1$ at $x_f$. If the condition $y_1(10)=200$ is not in the interval obtained by the two shootings the program stops.
4.2 - Subroutine to load initial vaues
Go to the topYou will notice that this subroutine is very small and having a whole subroutine only for this could be a waste of time and energy but for ODEs of higher order with larger number of conditions it could be useful. This is the code,
subroutine Load(yi)
use nrtypes
implicit none
real, intent(out):: yi(1:2)
!
! In this routine you may store all known boundary conditions at xi
yi(1) = 40.0_sp
!
return
end subroutine Load
4.3 - Subroutine for 4th-order Runge-Kutta method equations
Go to the topThis routine is in fact the same use before for solving initial value problems. Here it is,
subroutine rk4odes20240526(x,y,n,h)
use nrtypes
implicit none
real(SP), intent(inout) :: x
real(SP), intent(in) :: h
real(SP), dimension(:), intent(inout) :: y
integer(I4B), intent(in) :: n
real(SP), dimension(1:n) :: ym, ye, slope, yt
integer(I4B) :: nc ! Dummy counter for n
real(SP), DIMENSION(1:4,1:n) :: k
!
do nc = 1,n
yt(nc) = y(nc)
end do
!
! For k1
call derivs(x,y,1,k)
do nc = 1,n
ym(nc) = y(nc) + k(1,nc)*h/2.0
end do
!
! For k2
call derivs(x+h/2.0,ym,2,k)
do nc = 1,n
ym(nc) = y(nc) + k(2,nc)*h/2.0
end do
!
! For k3
call derivs(x+h/2.0,ym,3,k)
do nc = 1,n
ye(nc) = y(nc) + k(3,nc)*h
end do
!
! For k4
call derivs(x+h,ye,4,k)
!
! The slope is calculated as follows
do nc = 1,n
slope(nc) = ( k(1,nc) + 2*( k(2,nc) + k(3,nc) ) + k(4,nc) )/6.0
y(nc) = yt(nc) + slope(nc)*h
end do
x = x + h
return
end subroutine
In this routine the differential equations are evaluated according to the initial conditions provided and to the new values of y's estimated in 4th-order Runge-Kutta method routine. This is the code
subroutine derivs(x,y,kc,k)
use nrtypes
implicit none
REAL(SP), INTENT(IN) :: x
integer(I4B):: kc
REAL(SP), DIMENSION(:), INTENT(IN) :: y
REAL(SP), DIMENSION(6) :: dy
real(SP), DIMENSION(:,:), INTENT(OUT) :: k
!
dy(1) = y(2)
dy(2) = -0.01_sp*(20.0_sp-y(1))
!
k(kc,1) = dy(1)
k(kc,2) = dy(2)
return
end subroutine derivs
The program was executed and the output along with the results presented in the textbook by Chapra and Canale are shown below. You may try other values for the unknown condition $y_2$ at $x_i$. By the way the shooting for $y_2$ that gives $y_1(10)=200$ is 12.69.
 |
| Fig. 2 Comparison of results between those obtained from our program and those reported in the book of Chapra and Canale. |
There are several ideas to improve this program and it is left to the reader to implemented or adapted to his/her own programs.
- The arguments in the subroutines can be order so that the variables/parameters in it could be ordered according with an important feauture of these: input, input/output and output.
- You may modify the program to work with linear extrapolation rather than interpolation. You may use the same example to validate your results which should be the same.
- A subroutine dedicated to the shooting process can be taken from the main program so that it can be easily identified. As you can see, this is the heart of the program.
This is the end of the post. I hope you find it useful.
Any question? Write in the comments and I shall try to help.
Perhaps the following post may be of interest for you
- 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
- Dealing with higher order ODE when using the 4th-order Runge-Kutta mehtod
=========
Ildebrando
No comments:
Post a Comment