[Fortran] 纯文本查看 复制代码
program main
implicit none
integer, parameter :: nvar = 2, nstep=500
real(kind=8), dimension(nvar) :: ystart=[0.0d0,0.0d0]
real(kind=8), dimension(nstep) :: xx
real(kind=8), dimension(nvar,nstep) :: y, ytrue
real(kind=8) :: x1=0.0d0, x2=10.0d0
integer :: i
call rkdub(Qa,xx, y, ystart, nvar, x1, x2, nstep, diff) !这一行出现：error #6404: This name does not have a type, and must have an explicit type. [QA]
print*,y(1,i), y(2,i)
end program
! Define the Ordinary Differential Equations(ODEs) in this fileubroutine diff(Qa, t, y, dydx)
implicit none
integer, parameter :: n = 2
real(kind=8) :: Qa, t
real(kind=8) :: y(n), dydx(n)
real(kind=8) :: Kr=1.4
real(kind=8) :: Pg=5000000
real(kind=8) :: Vg=0.04
!real(kind=8), parameter :: Qa=0.0005
real(kind=8) :: Dm=0.00025
real(kind=8) :: Jm=2
real(kind=8) :: Ro=850
real(kind=8) :: Cd=0.61
real(kind=8) :: Af=0.0005
real(kind=8) :: Bg=2
! Define the ODE or ODEs here
!dx(1) = Kr*Pg*(Qa-Dm*x(2))/Vg
!dx(2) = (x(2)*x(2)*Ro*Dm*Dm*Dm/(2*(Cd*Af)*(Cd*Af))-Bg*x(2)+x(1)*Dm)/J
dydx(1) = Qa-Dm*y(2)
!dx(2) = (Dm*Pg/((1-x(1)/Vg)**1.4)-x(2)*x(2)*Ro*Dm*Dm*Dm/(2*(Cd*Af)*(Cd*Af))-Bg*x(2))/Jm
dydx(2) = (Dm*(Vg**Kr*Pg/((Vg-y(1))**1.4))-(y(2)*y(2)*Ro*Dm*Dm*Dm)/(2*(Cd*Af)*(Cd*Af))-Bg*y(2))/Jm
end subroutine diff
! RK4 numerical solution of differential equations
subroutine ODE_RK4(Qa,t, y, n, h, yout, diff)
implicit none
external :: diff
integer n
real(kind=8) :: Qa, t, h
real(kind=8) :: y(n)
real(kind=8) :: yout(n)
real(kind=8) :: k1(n),k2(n),k3(n),k4(n)
! Calculate k1
call diff(Qa, t, y, dydx) ！这一行出现：error #6404: This name does not have a type, and must have an explicit type. [QA]
k1 = h*dydx
call diff(Qa, t+0.5d0*h, y+0.5*k1, dydx)
k2 = h*dydx
call diff(Qa, t+0.5d0*h, y+0.5*k2, dydx)
k3 = h*dydx
call diff(Qa, t+h, y+k3, dydx)
k4 = h*dydx
yout = y + k1/6.0d0 + k2/3.0d0 + k3/3.0d0 + k4/6.0d0
end subroutine ODE_RK4
subroutine rkdub(Qa,xx, y, ystart, nvar, x1, x2, nstep, diff)
implicit none
integer, intent(in) :: nvar, nstep, Qa
real(kind=8), intent(in) :: x1, x2
real(kind=8), dimension(nvar), intent(in) :: ystart
real(kind=8), dimension(nstep), intent(out) :: xx
real(kind=8),dimension(nvar,nstep), intent(out) :: y
external :: diff
real(kind=8) :: h, x
integer :: i,j
x = x1
h = (x2-x1) / nstep
xx(1) = x1
y(:,1) = ystart
do i=2, nstep
call ODE_RK4(Qa,x, y(:,i-1), nvar, h, y(:,i), diff)
if (x+h == x) then
write(*,*) 'stepsize not significant in rkdub'
end if
x = x + h
xx(i) = x
end do
end subroutine rkdub