[Fortran] 纯文本查看 复制代码
module counts
integer,parameter::nstep=100
integer,parameter::nvar=3 !Funcion dimentions
end module counts
subroutine derivs(n,x,y,dydx)
implicit none
integer::n !方程维数
real::dydx(N),x(N),y(N)
dydx(1)=y(2)*y(3)
dydx(2)=-y(1)*y(3)
dydx(3)=-0.51*y(1)*y(2)
end subroutine derivs
subroutine rk4(y,dydx,n,x,h,yout,derivs)
integer,parameter::nmax=3
real::h,x,dydx(n),y(n),yout(n)
external::derivs
integer::i
real::h6,hh,xh,dym(nmax),dyt(nmax),yt(nmax)
hh=h*0.5
h6=h/6.
xh=x+hh
do i=1,n !First step,k1=dydx(i)
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs(xh,yt,dyt) !Second step,k2=dyt
do i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs(xh,yt,dym) !Third step,k3=dym
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i) !k3=k3+k2
enddo
call derivs(x+h,yt,dyt)!Fourth step,k4=dyt
do i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
end do !Every step calculate dydx from x and above y(i)
return
end subroutine rk4
subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs)
integer::nstep,nvar
integer,parameter::nmax=3,nstpmx=100
real::x2,x1,vstart(nvar),xx(nstpmx),y(nmax,nstpmx)
external derivs
integer::i,k
real::h,x,dv(nmax),v(nmax)
do i=1,nvar !方程维数
v(i)=vstart(i)
y(i,1)=v(i)
enddo
xx(1)=x1
x=x1
h=(x2-x1)/nstep !步长
do k=1,nstep
call derivs(x,v,dv)
call rk4(v,dv,nvar,x,h,v,derivs)
if(x+h.eq.x) pause'stepsize not significant in rkdumb.'
x=x+h
xx(k+1)=x
do i=1,nvar
y(i,k+1)=v(i)
enddo
enddo
return
end subroutine rkdumb
program solve319
use counts
implicit none
integer::i
integer,parameter::nmax=3,nstpmx=100
real::x2,x1,vstart(nvar),xx(nstpmx),y(nmax,nstpmx)
external::derivs,rk4,rkdumb
vstart(1)=0
vstart(2)=1
vstart(3)=1
x1=0
x2=12
print*,vstart(1)
open(unit=11,file='20150327out.dat')
do i=1,nstep
write(unit=11,fmt='(3f12.7)')y(1,i),y(2,i),y(3,i)
call rkdumb(vstart,nvar,x1,x2,nstep,derivs)
enddo
end