【2階微分　ルンゲクッタ法にて】


csvファイル出力をしない版はこちら
****** main ******
      external f1,f2
      
      write(*,*)'h='
      read(*,*)h
      
      write(*,*)'t0='
      read(*,*)t0
      
      write(*,*)'y0='
      read(*,*)y0
      
      write(*,*)'y0''='
      read(*,*)y0d
      
      write(*,*)'tend='
      read(*,*)tend
      
      call rk2(f1,f2,h,t0,y0,y0d,tend)
      stop
      end
      
****** subroutine program*******
      subroutine rk2(f1,f2,h,t0,y0,y0d,tend)
      real*4 j1,j2,j3,j4,k1,k2,k3,k4
      t =t0
      y =y0
      yd=y0d
      open(7,file='runge2.csv', status='replace')
      write(7,*)'h,',h
      write(7,*)'t,y,yd,j1,j2,j3,j4,k1,k2,k3,k4'
      
      write(7,*)t,',',y,',',yd
   20 continue
      if((t+0.0001).gt.tend)go to 30
      j1= f1(y,yd,t)
      k1= f2(y,yd,t)
      j2= f1(y+h*j1/2,yd+h*k1/2,t+h/2)
      k2= f2(y+h*j1/2,yd+h*k1/2,t+h/2)
      j3= f1(y+h*j2/2,yd+h*k2/2,t+h/2)
      k3= f2(y+h*j2/2,yd+h*k2/2,t+h/2)
      j4= f1(y+h*j3,yd+h*k3,t+h)
      k4= f2(y+h*j3,yd+h*k3,t+h)
      y = y+(h/6)*(j1+2*j2+2*j3+j4)
      yd=yd+(h/6)*(k1+2*k2+2*k3+k4)
      t =t+h
      write(7,*)t,',',y,',',yd,',',j1,',',j2,',',j3,',',j4
     a,',',k1,',',k2,',',k3,',',k4
      goto 20
      
   30 write(*,*)' y  =',y
      write(7,*)'y,',y
      close(7)
      return
      end
      
****** function f1 *******
      function f1(y,yd,t)
      f1 = yd
      return
      end
      
****** function f2 *******
      function f2(y,yd,t)
      f2 = t + 1
      return
      end