29 lines
926 B
Forth
29 lines
926 B
Forth
program rungekutta
|
|
implicit none
|
|
real(kind=kind(1.0D0)) :: t,dt,tstart,tstop
|
|
real(kind=kind(1.0D0)) :: y,k1,k2,k3,k4
|
|
tstart =0.0D0 ; tstop =10.0D0 ; dt = 0.1D0
|
|
y = 1.0D0
|
|
t = tstart
|
|
write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&
|
|
&,abs(y-(t**2+4.0d0)**2/16.0d0)
|
|
do; if ( t .ge. tstop ) exit
|
|
k1 = f (t , y )
|
|
k2 = f (t+0.5D0 * dt, y +0.5D0 * dt * k1)
|
|
k3 = f (t+0.5D0 * dt, y +0.5D0 * dt * k2)
|
|
k4 = f (t+ dt, y + dt * k3)
|
|
y = y + dt *( k1 + 2.0D0 *( k2 + k3 ) + k4 )/6.0D0
|
|
t = t + dt
|
|
if(abs(real(nint(t))-t) .le. 1.0D-12) then
|
|
write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&
|
|
&,abs(y-(t**2+4.0d0)**2/16.0d0)
|
|
end if
|
|
end do
|
|
contains
|
|
function f (t,y)
|
|
implicit none
|
|
real(kind=kind(1.0D0)),intent(in) :: y,t
|
|
real(kind=kind(1.0D0)) :: f
|
|
f = t*sqrt(y)
|
|
end function f
|
|
end program rungekutta
|