CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C THE FOLLOWING SUBROUTINES ARE FOR FIXED STEPSIZE C EXPLICIT RUNGE-KUTTA METHODS OF ORDERS FROM 1 TO 6. C THESE SUBROUTINES ARE TO CALCULATE NUMERICAL SOLUTION OF AN C AUTONOMOUS SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL C EQAUTIONS Y'=F(Y). C CC VARIABLE AND ARRAY: C YN: INPUT AND OUTPUT ARRAY OF VALUES FOR APPROXIMATE SOLUTION; C STORE INITIAL VALUES AND OUTPUT APPROXIMATE VALUES AFTER C ONE INTEGRATION STEP. C N: DIMENSION OF THE SYSTEM ( I.E., NUMBER OF EQUATIONS, C OR NUMBER OF VARIABLES) C H: INTEGRATION STEPSIZE C C SUBROUTINES: C subroutine eu1(n,yn,h) : the Euler's method C subroutine rk2(n,yn,h) : the improved Euler's method C subroutine rk3(n,yn,h) : a Runge-Kutta method of order 3 C subroutine rk4(n,yn,h) : a Runge-Kutta method of order 4 C subroutine rk5(n,yn,h) : a Runge-Kutta method of order 5 C subroutine rk6(n,yn,h) : a Runge-Kutta method of order 6 C subroutine rkm2(n,yn,h): another Runge-Kutta method of order 2 C subroutine rkh3(n,yn,h): another Runge-Kutta method of order 3 C CC SUBROUTINE DIFFUN COMPILED BY YOURSELF: C DIFFUN IS NAME OF SUBROUTINE COMPUTING THE FIRST DERITAVE F(Y), C AND THIS SUBROUTINE IS OF THE FORM C SUBROUTINE DIFFUN(N,Y,DY) C DIMENSION Y(N),DY(N) C Y: INITIAL VALUES C N: NUMBER OF VARIABLES C DY: FIRST DERITAVES F(Y) C FOR EXAMPLE: C FOR THE FIRST ORDER ORDINARY DIFFERENTIAL EQAUTIONS C Y1'=-Y2 C Y2'=Y1 C C subroutine diffun(n,y,dy) C dimension y(n),dy(n) C dy(1)=-y(2) C dy(2)=y(1) C return C end C COLLECTED BY Dr. Jianping Li, DEC 20, 1998. c------------------------------------------------------------ c This subroutine is for the Euler's method c which is a Runge-Kutta method of order one. c By Dr. Jianping Li, Sep. 1, 1997. c----- subroutine eu1(n,yn,h) dimension yn(n) dimension dy(n) !work array call diffun(n,yn,dy) do 10 i=1,n 10 yn(i)=yn(i)+h*dy(i) return end c------------------------------------------------------------ c This subroutine is for the improved Euler's method c which is a Runge-Kutta method of order 2. c By Dr. Jianping Li, Sep. 3, 1997. c----- subroutine rk2(n,yn,h) parameter(m=2) dimension yn(n) dimension y0(n),y(n),dk(n,m) !work array do 10 i=1,n y0(i)=yn(i) 10 continue call diffun(n,yn,dk(1,1)) do 20 i=1,n y(i)=y0(i)+h*dk(i,1) 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n ay=(dk(i,1)+dk(i,2))/2. yn(i)=yn(i)+h*ay 30 continue return end c------------------------------------------------------------ c This subroutine is for the Kutta method c which is a Runge-Kutta method of order 3. c By Dr. Jianping Li, Sep. 5, 1997. c----- subroutine rk3(n,yn,h) parameter(m=3) dimension yn(n) dimension u(2),y0(n),y(n),dk(n,m) !work array u(1)=h/2.0 u(2)=2.*h do 10 i=1,n y0(i)=yn(i) y(i)=yn(i) 10 continue call diffun(n,y,dk(1,1)) do 20 i=1,n y(i)=y0(i)+u(1)*dk(i,1) 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n y(i)=y0(i)-h*dk(i,1)+u(2)*dk(i,2) 30 continue call diffun(n,y,dk(1,3)) do 40 i=1,n ay=(dk(i,1)+4.*dk(i,2)+dk(i,3))/6. yn(i)=yn(i)+h*ay 40 continue return end c------------------------------------------------------------ c This subroutine is for the classical c Runge-Kutta method of order 4. c By Dr. Jianping Li, Sep. 5, 1997. c----- subroutine rk4(n,yn,h) parameter(m=4) dimension yn(n) dimension u(3),y0(n),y(n),dk(n,m) !work array u(1)=h/2.0 u(2)=u(1) u(3)=h do 10 i=1,n y0(i)=yn(i) y(i)=yn(i) 10 continue do 30 j=1,3 call diffun(n,y,dk(1,j)) do 20 i=1,n y(i)=y0(i)+u(j)*dk(i,j) 20 continue 30 continue call diffun(n,y,dk(1,4)) do 40 i=1,n ay=(dk(i,1)+2.*(dk(i,2)+dk(i,3))+dk(i,4))/6. yn(i)=yn(i)+h*ay 40 continue return end c------------------------------------------------------------ c This subroutine is for the Fehlberg's 5-th order method c which is a 6-stage 5-th order Runge-Kutta method. c By Dr. Jianping Li, Oct. 15, 1997. c----- subroutine rk5(n,yn,h) parameter(m=6) dimension yn(n) dimension y0(n),y(n),dk(n,m) !work array do 10 i=1,n y0(i)=yn(i) y(i)=yn(i) 10 continue call diffun(n,y,dk(1,1)) do 20 i=1,n y(i)=y0(i)+h*dk(i,1)/4. 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n y(i)=y0(i)+h*(3.*dk(i,1)+9.*dk(i,2))/32. 30 continue call diffun(n,y,dk(1,3)) do 40 i=1,n y(i)=y0(i)+h*(1932.*dk(i,1)-7200.*dk(i,2) * +7296.*dk(i,3))/2197. 40 continue call diffun(n,y,dk(1,4)) do 50 i=1,n y(i)=y0(i)+h*(8341.*dk(i,1)-32832.*dk(i,2) * +29440.*dk(i,3)-845.*dk(i,4))/4104. 50 continue call diffun(n,y,dk(1,5)) do 60 i=1,n y(i)=y0(i)+h*((-1216.*dk(i,1)+1859.*dk(i,4))/4104. * +2.*dk(i,2)-3544.*dk(i,3)/2565. * -11.*dk(i,5)/40.) 60 continue call diffun(n,y,dk(1,6)) do 70 i=1,n ay=16.*dk(i,1)/135.+6656.*dk(i,3)/12825. * +28561.*dk(i,4)/56430. * -9.*dk(i,5)/50.+2.*dk(i,6)/55. yn(i)=yn(i)+h*ay 70 continue return end c------------------------------------------------------------ c This subroutine is for the Butcher's 6-th order method c which is a 6-th order Runge-Kutta method with 7 stages. c By Dr. Jianping Li, Oct. 18, 1997. c----- subroutine rk6(n,yn,h) parameter(m=7) dimension yn(n) dimension y0(n),y(n),dk(n,m) !work array do 10 i=1,n y0(i)=yn(i) y(i)=yn(i) 10 continue call diffun(n,y,dk(1,1)) do 20 i=1,n y(i)=y0(i)+h*dk(i,1)/2. 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n y(i)=y0(i)+h*(2.*dk(i,1)+4.*dk(i,2))/9. 30 continue call diffun(n,y,dk(1,3)) do 40 i=1,n y(i)=y0(i)+h*(7.*dk(i,1)+8.*dk(i,2)-3.*dk(i,3))/36. 40 continue call diffun(n,y,dk(1,4)) do 50 i=1,n y(i)=y0(i)+h*(-35.*dk(i,1)-220.*dk(i,2)+105.*dk(i,3) * +270.*dk(i,4))/144. 50 continue call diffun(n,y,dk(1,5)) do 60 i=1,n y(i)=y0(i)+h*(-dk(i,1)-110.*dk(i,2)-45.*dk(i,3) * +180.*dk(i,4)+36.*dk(i,5))/360. 60 continue call diffun(n,y,dk(1,6)) do 70 i=1,n y(i)=y0(i)+h*(-41.*dk(i,1)/260.+22.*dk(i,2)/13. * +43.*dk(i,3)/156.-118.*dk(i,4)/39. * +32.*dk(i,5)/195.+80.*dk(i,6)/39.) 70 continue call diffun(n,y,dk(1,7)) do 120 i=1,n ay=(13.*(dk(i,1)+dk(i,7))+55.*(dk(i,3)+dk(i,4)) * +32.*(dk(i,5)+dk(i,6)))/200. yn(i)=yn(i)+h*ay 120 continue return end c------------------------------------------------------------ c This subroutine is for the midpoint method c which is a Runge-Kutta method of order 2. c By Dr. Jianping Li, Nov. 1, 1997. c----- subroutine rkm2(n,yn,h) parameter(m=2) dimension yn(n) dimension y0(n),y(n),dk(n,m) !work array do 10 i=1,n y0(i)=yn(i) 10 continue call diffun(n,yn,dk(1,1)) do 20 i=1,n y(i)=y0(i)+h*dk(i,1)/2. 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n yn(i)=yn(i)+h*dk(i,2) 30 continue return end c------------------------------------------------------------ c This subroutine is for the Heun method c which is a Runge-Kutta method of order 3. c By Dr. Jianping Li, Nov. 3, 1997. c----- subroutine rkh3(n,yn,h) parameter(m=3) dimension yn(n) dimension u(2),y0(n),y(n),dk(n,m) !work array do 10 i=1,n y0(i)=yn(i) y(i)=yn(i) 10 continue call diffun(n,y,dk(1,1)) do 20 i=1,n y(i)=y0(i)+h*dk(i,1)/3. 20 continue call diffun(n,y,dk(1,2)) do 30 i=1,n y(i)=y0(i)+2.*h*dk(i,2)/3. 30 continue call diffun(n,y,dk(1,3)) do 40 i=1,n ay=(dk(i,1)+3.*dk(i,3))/4. yn(i)=yn(i)+h*ay 40 continue return end