简单来说就是解微分方程
用4阶的Runge Kutta
程式码如下
SUBROUTINE RKINT (DERIVY, Y, DY, Q, NEQFST, NEQLST, DX)
C ref. Ralston & Wilf "Mathematical Methods for Digital Computers",p.117
REAL*8 Y(NEQLST), DY(NEQLST), Q(NEQLST)
INTEGER NEQFST,NEQLST
REAL*8 A(4), B(4), C(4)
REAL*8 DX,T
EXTERNAL DERIVY
DATA A /0.5D0, 0.29289322D0, 1.7071068D0, 0.16666667D0/,
1 B /2.0D0, 1.0D0, 1.0D0, 2.0D0/,
2 C /0.5D0, 0.29289322D0, 1.7071068D0, 0.5D0/
DO 1 J = 1, 4
CALL DERIVY (Y, DY, NEQFST, NEQLST)
DO 2 I = NEQFST, NEQLST
T = A(J)*(DY(I) - B(J)*Q(I))
Y(I) = Y(I) + DX*T
Q(I) = Q(I) + 3.0D0*T - C(J)*DY(I)
2 CONTINUE
1 CONTINUE
RETURN
END
其中DERIVY是描述微分方程的副程式,NEQLST是方程式的数目
Y是变量值 DY是变量的微分值
目前解的问题有8个变量8个方程式
但是今天发现如果DY太大会出现array bound exceed的情况...
array bound exceed我知道是使用的指令超出矩阵宣告范围
我在DERIVY中是将DY(1)到DY(8)都设10000，就出现array bound exceed
设小一点都没事 实在想不透原因...拜托版上强者帮忙了...