簡單來說就是解微分方程
用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
設小一點都沒事 實在想不透原因...拜託版上強者幫忙了...