Quellcode-Bibliothek odeint.ftn Sprache: Fortran

c<html>
c<head><title>odeint.f</title></head>
c<body>
c<pre>
c<a name="module">
 module precision
c</a>
c<a name=2>
 integer, parameter :: rp=selected_real_kind(13,300)
c
 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">
 1 write(*,'(a)',advance='no') 'Provide Integration Step Size: '
c </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>
 if(abs(dely/max(yimp,yimpn)).lt.1.e-9) go to 72
c ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 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">
 end function func
c</a>
 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>
 e2x=exp(2*x)
c
 y=(e2x-1.)/(1.+e2x)
 return
 end
c</pre>
c</body>
c</html>

quality79%

¤ Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.0.13Bemerkung: (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.