!************************************************************************** !* Solve a first order DE system (N=2) of the form: * !* y' = F(x,y,z), z'=G(x,y,z) with conditions: y(x0)=a and z(x0)=b from * !* x=x0 to x=x1 using a Runge-Kutta integration method. * !* ---------------------------------------------------------------------- * !* SAMPLE RUN: * !* Integrate first order DE systemy: * !* y' = y*z + cos(x) - 0.5*sin(2*x) * !* z' = y*y + z*z - (1 + sin(x)) * !* from x=0 to x=2 with initial conditions: y(0)=0 and z(0)=0 * !* (25 integration steps). * !* * !* X Y estimated Z estimated * !* ------------------------------------ * !* 0.0000 0.0000000 0.0000000 * !* 0.0800 0.0765518 -0.0828580 * !* 0.1600 0.1452904 -0.1700792 * !* 0.2400 0.2050556 -0.2597237 * !* 0.3200 0.2550329 -0.3500533 * !* 0.4000 0.2948004 -0.4396008 * !* 0.4800 0.3243379 -0.5271915 * !* 0.5600 0.3440043 -0.6119278 * !* 0.6400 0.3544874 -0.6931461 * !* 0.7200 0.3567373 -0.7703644 * !* 0.8000 0.3518908 -0.8432315 * !* 0.8800 0.3411943 -0.9114874 * !* 0.9600 0.3259299 -0.9749377 * !* 1.0400 0.3073506 -1.0334409 * !* 1.1200 0.2866242 -1.0869048 * !* 1.2000 0.2647908 -1.1352886 * !* 1.2800 0.2427308 -1.1786041 * !* 1.3600 0.2211465 -1.2169165 * !* 1.4400 0.2005537 -1.2503410 * !* 1.5200 0.1812830 -1.2790362 * !* 1.6000 0.1634888 -1.3031952 * !* 1.6800 0.1471644 -1.3230351 * !* 1.7600 0.1321612 -1.3387869 * !* 1.8400 0.1182104 -1.3506867 * !* 1.9200 0.1049465 -1.3589674 * !* 2.0000 0.0919314 -1.3638530 * !* * !* ---------------------------------------------------------------------- * !* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en Basic * !* et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"* !* [BIBLI 04]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !************************************************************************** Program Test_Rk4 integer, parameter :: SIZE = 100 real*8 X(SIZE),Y(SIZE),Z(SIZE) real*8 xl,x0,h integer k,kl,l pause ' Adjust window size...' x0=0.d0 !starting x xl=2.d0 !ending x kl=25 !number of steps in x h=(xl-x0)/kl !integration step X(0)=x0 Y(0)=0.d0; Z(0)=0.d0 !starting values !integration loop do k=0, kl-1 call RK4(X(k),Y(k),Z(k),h,X(k+1),Y(k+1),Z(k+1)) end do !write header print *,' ' print *,' X Y estimated Z estimated ' print *,' -------------------------------------' !write kl+1 result lines do k=0, kl write(*,10) X(k), Y(k), Z(k) end do print *,' -------------------------------------' print *,' ' stop 10 format(F9.4,F12.7,F15.7) END !y'=yz + cos(x) - 0.5sin(2x) real*8 Function F(x,y,z) real*8 x,y,z F=y*z + dcos(x) - 0.5d0*dsin(2.d0*x) return End !z'=yy + zz -(1+sin(x)) real*8 Function G(x,y,z) real*8 x,y,z G=y*y + z*z -(1.d0+dsin(x)) return End !Integrate sytem from x to x+h using Runge-Kutta Subroutine RK4(x,y,z,h,x1,y1,z1) real*8 x,y,z,h,x1,y1,z1 real*8 c1,c2,c3,c4,d1,d2,d3,d4,h2,F,G c1=F(x,y,z) d1=G(x,y,z) h2=h/2.d0 c2=F(x+h2,y+h2*c1,z+h2*d1) d2=G(x+h2,y+h2*c1,z+h2*d1) c3=F(x+h2,y+h2*c2,z+h2*d2) d3=G(x+h2,y+h2*c2,z+h2*d2) c4=F(x+h,y+h*c3,z+h*d3) d4=G(x+h,y+h*c3,z+h*d3) x1=x+h y1=y+h*(c1+2.d0*c2+2.d0*c3+c4)/6.d0 z1=z+h*(d1+2.d0*d2+2.d0*d3+d4)/6.d0 return End !end of file trk4.f90