c<html>
c<head><title>odeint.f</title></head>
c<body>
c<pre>
c<a name="module"><font color="FF0900">
module precision
c</font></a>
c<a name=2><font color=FF0000>
integer, parameter :: rp=selected_real_kind(13,300)
c</font>
end module precision
program odeint
use precision
implicit none
c
c Demonstration of Runge-Kutta, Adams-Bashford, and Adams-Moulton for
c solution of an Ordinary Differential Equation in the form:
c
c d y
c --- = f ( x , y )
c d x
c
c The arrays below allow storage of the results for the four most recent
c time steps. For example y(0) contains the value of y at the beginning
c of the current time step, y(1) contains the value of y at the beginning
c or the previous time step, y(2) contains the value of y 2 steps back, etc.
c
real (rp) f(0:3),y(0:3),fa(0:3),ya(0:3),fimp(0:3),yt(0:3),ft(0:3)
real (rp) h, ynew, yrk, yab, yanew, xnew, dfdy, xn, xk1,
& xk2, xk3, xk4, yam, yimpn, yimp, gy, dgdy, dely, yabt,
& fnew, func, ftnew
integer nsteps, istep, i, it
character*12 filename
character*5 xstep
c Read in the integration step size
c <a name="write"><font color="FF0000">
1 write(*,'(a)',advance='no') 'Provide Integration Step Size: '
c </font></a>
read (*,*)h
c Make up output a file name that means something
c Combining 'odeout-' with the step size
write(xstep,'(f5.3)') h
filename = 'odeout-'//adjustl(xstep)
open(10,file=filename)
write(10,2010)
2010 format(4x,'x',10x,'correct',5x,'Runge',7x,'Adams',7x,'Adams',4x,
$ 'Implicit',/15x,'answer',6x,'Kutta',7x,'Bashford',4x,'Moulton',
$ 2x,'Adams-Moultn',/)
c
c Get A number of integration steps
c
write(*,'(a)',advance='no') 'Number of integration steps: '
read(*,*) nsteps
c
c Take at least 5 steps
c
nsteps=max(nsteps,5)
c
c Here We hardwire the initial conditions
c
ynew=0.
xnew=0.
call answer(xnew,yanew)
c
c Initialize values of previous timesteps
c
do 10 i=1,3
y(i)=0.
ya(i)=0.
yt(i)=0.
f(i)=0.
fimp(i)=0.
ft(i)=0.
10 fa(i)=0.
c
c Do Three steps of Runge Kutta
c
do 40 istep=1,3
c
c Shuffle past data
c
do 30 i=3,1,-1
yt(i)=yt(i-1)
ya(i)=ya(i-1)
y(i)=y(i-1)
f(i)=f(i-1)
ft(i)=ft(i-1)
fimp(i)=fimp(i-1)
30 continue
y(0)=ynew
yt(0)=ynew
ya(0)=yanew
xn= xnew
c
c Make evaluations for this time step
c
f(0)=func(xn,y(0),dfdy)
ft(0)=f(0)
c
c Initialize numbers needed for the Implicit Adams-Moulton method.
c Note that I am cheating here fimp should be
c loaded with the results of an implicit Runge-Kutta
c or explicit formulation for variable step size using
c smaller steps for initialization
c
fimp(0)=func(xn,ya(0),dfdy)
c
c Actual Runge-Kutta step
c
xk1=h*func(xn,y(0),dfdy)
xk2=h*func(xn+.5*h,y(0)+.5*xk1,dfdy)
xk3=h*func(xn+.5*h,y(0)+.5*xk2,dfdy)
xk4=h*func(xn+h,y(0)+xk3,dfdy)
c
c Values at the new time level (end of time step)
c If you look at this carefully it is a close relative of a
c Simpson's Rule integration
c
ynew=y(0)+(xk1+2.*xk2+2.*xk3+xk4)/6.
xnew=xn+h
c
c Get the actual answer for comparison
c
call answer(xnew,yanew)
write(10,2000)xnew,yanew,ynew
40 continue
2000 format(1p,6e12.5)
c
c Now that we have starter values for Adams-Bashford and
c Adams-Moulton, Finish integration with all methods
c
c Current y for Runge-Kutta
c
yrk=ynew
c
c Current y for Adams-Bashford
c
yab=ynew
c
c Current y for Adams-Moulton
c
yam=ynew
c
c Current y for Implicit Adams-Moulton
c
yimpn=yanew
c
c
c
do 80 istep=4,nsteps
c
c Set up for next time step by shifting stored values of previous time
c steps to keep only the results of the three previous steps
c
do 60 i=3,1,-1
ya(i)=ya(i-1)
yt(i)=yt(i-1)
y(i)=y(i-1)
f(i)=f(i-1)
ft(i)=ft(i-1)
fimp(i)=fimp(i-1)
60 continue
xn= xnew
y(0)=yam
yt(0)=yab
f(0)=func(xn,y(0),dfdy)
ft(0)=func(xn,yt(0),dfdy)
c
c The Implicit Adams-Moulton requires a Newton Iteration (DO loop 70)
c to solve a non-linear equation in the new-time value of y
c
yimp=yimpn
fimp(0)=func(xn,yimp,dfdy)
xnew=xn+h
do 70 it=1,10
fnew=func(xnew,yimpn,dfdy)
gy=yimp-yimpn+(9.*fnew+19.*fimp(0)-5.*fimp(1)+fimp(2))*h/24.
dgdy=9.*h/24.*dfdy-1.
dely=-gy/dgdy
yimpn=yimpn+dely
c <a name=abs><font color="FF0000">
if(abs(dely/max(yimp,yimpn)).lt.1.e-9) go to 72
c </font> ^^^^^^^^^^^^^^^^^^^^^^^^^^^
70 continue
write(6,*) 'Implicit iteration failed'
stop
c
c Standard Runge-Kutta continues
c
72 xk1=h*func(xn,yrk,dfdy)
xk2=h*func(xn+.5*h,yrk +.5*xk1,dfdy)
xk3=h*func(xn+.5*h,yrk +.5*xk2,dfdy)
xk4=h*func(xn+h,yrk +xk3,dfdy)
yrk =yrk +(xk1+2.*xk2+2.*xk3+xk4)/6.
c
c Do a straight Adams-Bashford integration
c
yab=yam+(55.*ft(0)-59.*ft(1)+37.*ft(2)-9.*ft(3))*h/24.
c
c Do an Adams-Bashford predictor for Adams-Moulton
c
yabt=yam+(55.*f(0)-59.*f(1)+37.*f(2)-9.*f(3))*h/24.
ftnew=func(xnew,yabt,dfdy)
c
c Now apply and Adams-Moulton correction step
c
yam=yam+(9.*ftnew+19.*f(0)-5.*f(1)+f(2))*h/24.
call answer(xnew,yanew)
write(10,2000)xnew,yanew,yrk,yab,yam,yimpn
80 continue
c
close(10)
stop
end
c
function func(x,y,dfdy)
use precision
real (rp) x,y,dfdy,func
c
c A specific function f for the right hand side of the ODE
c
func=1.-y**2
c
c Limit the value to keep the code running when one numerical solution
c method goes nuts
c
func=min(1.e10_rp,max(func,-1.e10_rp))
dfdy=-2*y
return
c<a name="ef"><font color="FF0000">
end function func
c</a></font>
subroutine answer(x,y)
use precision
real (rp) x,y,e2x
c
c For the function in "func" this is the analytic solution to the ODE
c<a name=3><font color=FF0000>
e2x=exp(2*x)
c</font>
y=(e2x-1.)/(1.+e2x)
return
end
c</pre>
c</body>
c</html>
Copyright ©2017 - STI - Todos os direitos reservados