!******************************************************************** !* Differential equations of order N * !* by Runge-Kutta method of order 4 * !* with boundary conditions * !* ---------------------------------------------------------------- * !* Reference: 'Analyse en Turbo Pascal versions 5.5 et 6.0' by Marc * !* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991. * !* * !* F90 Version By J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* ---------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6 * !* with initial conditions: y(4)=2 and y'(4)=-2 * !* and final condition y(6)=2 * !* * !* Exact solution is: y = 2cos(1)/cos(x-5) * !* * !* DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N * !* of type y(n) = f(y(n-1),...,y',y,x) * !* with boundary conditions * !* * !* Order: 2 * !* Begin value x : 4.00000000000000 * !* end value x : 6.00000000000000 * !* Y value at begin x : 2.00000000000000 * !* Y' value at begin x : -2.00000000000000 * !* Desired value at end x: 2.00000000000000 * !* number of points : 11 * !* finesse : 20 * !* Max. iterations : 100 * !* * !* Last iteration: * !* u =-3.11481552124023 du =-7.629394531250000E-007 * !* dyb= 8.195419964884820E-008 67 steps * !* * !* X Y * !* --------------------------- * !* 4.000000 2.000000 * !* 4.200000 1.551018 * !* 4.400000 1.309291 * !* 4.600000 1.173217 * !* 4.800000 1.102583 * !* 5.000000 1.080605 * !* 5.200000 1.102583 * !* 5.400000 1.173217 * !* 5.600000 1.309291 * !* 5.800000 1.551018 * !* 6.000000 2.000000 * !* * !******************************************************************** PROGRAM TSHOOT integer fi real*8 vect(1:50), xi,xf,yi(10) real*8 du,dyb,olddyb,u,yb print *,' ' print *,' DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N' print *,' of type y(n) = f(y(n-1),...,y'',y,x)' print *,' with boundary conditions' print *,' ' iordre=2 print *,' Order: ', iordre xi=4.d0; xf=6.d0 print *,' Begin value x: ', xi print *,' End value x : ', xf yi(0)=2.d0; yi(1)=-2.d0 print *,' Y value at begin x :', yi(0) print *,' Y'' value at begin x :', yi(1) yb=2.d0 print *,' Desired value at end x:', yb ndata=11 print *,' Number of points:', ndata fi=20 print *,' Finesse :', fi itermax=100 print *,' Max. iterations :', itermax vect=0.d0 u=yi(1) !First derivative initial value du=dabs(u)/10.d0 niter=1 ! estimate initial error on desired final value call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.false.) olddyb=dabs(vect(ndata)-yb) dyb=olddyb ! Iterate until desired final value yb is obtained do while(dyb>1.d-7.and.niter<=itermax) u=u+du; yi(1)=u call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.false.) dyb=dabs(vect(ndata)-yb) if (dabs(dyb)>dabs(olddyb)) then !overshoot? if (niter==1) then du=-du else du=-du/2.d0 end if end if olddyb=dyb niter=niter+1 end do print *,' ' print *,' Last iteration:' print *,' u =',u,' du =',du print *,' dyb=',dyb,' ',niter,' steps' !call subroutine equadifn() and print results call Equadifn(vect,xi,xf,yi,ndata,iordre,fi,.true.) stop 50 format(' y',i1,' value at x0 : ') END !Example: y'=(2y'y'+yy)/y real*8 FUNCTION fp(x,y) real*8 x,y(10) if (dabs(y(0))<1.d-12) y(0)=1.d-12 fp=(2.d0*y(1)*y(1)+y(0)*y(0))/y(0) end !*************************************************************************** !* SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N * !* of type y(n) = f(y(n-1),y(n-2),...,y',y,x) * !* ----------------------------------------------------------------------- * !* INPUTS: * !* xi, xf Begin, end values of variable x * !* yi Begin values at xi (of f and derivatives) * !* m Number of points to calculate * !* n Order of differential equation * !* fi finesse (number of intermediary points) * !* * !* OUTPUTS: * !* t real vector storing m results for function y * !* ----------------------------------------------------------------------- * !* EXAMPLE: y' = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1) * !* Exact solution: y = (2 cos(1))/cos(x-5) * !*************************************************************************** Subroutine Equadifn(t,xi,xf,yi,m,n,fi,flag) integer fi logical flag real*8 t(*),xi,xf,yi(*), fp real*8 h,x,a,b,c,d real*8 ta(10),tb(10),tc(10),td(10),y(10),z(10) ! print *,' ** m=', m ! print *,' ** ordre=', n ! print *,' ** fi=', fi if (fi<1) return h = (xf - xi) / fi / dfloat(m-1) n1=n-1; n2=n-2 t(1)=yi(0) do k=0, n1 y(k)=yi(k); z(k)=yi(k) end do do i=1, m ni=(i-1)*fi-1 do j=1, fi x=xi+h*(ni+j) do k=0, n1 y(k)=z(k) end do a=h*fp(x,y) do k=0, n2 ta(k)=h*y(k+1); y(k)=z(k)+ta(k)/2.d0 end do y(n1)=z(n1)+a/2.d0 x=x+h/2 b=h*fp(x,y) do k=0, n2 tb(k)=h*y(k+1); y(k)=z(k)+tb(k)/2.d0 end do y(n1)=z(n1)+b/2.d0 c=h*fp(x,y) do k=0, n2 tc(k)=h*y(k+1); y(k)=z(k)+tc(k) end do y(n1)=z(n1)+c x=x+h/2 d=h*fp(x,y) do k=0, n2 z(k)=z(k)+(ta(k)+2.d0*tb(k)+2.d0*tc(k)+h*y(k+1))/6.d0 end do z(n1)=z(n1)+(a+b+b+c+c+d)/6.d0 end do t(i+1)=z(0) end do !call subroutine Display if required if (flag) call Display(m,xi,xf,t) return end !print solutions Subroutine Display(m,xi,xf,t) integer i,m real*8 h,x,xi,xf,t(*) h=(xf-xi)/dfloat(m-1) x=xi-h print *,' ' print *, ' X Y ' print *, '----------------------' do i=1, m x=x+h write(*,100) x, t(i) end do print *,' ' return 100 format(' ',f9.6,' ',f9.6) end !End of file tshoot.f90