作者deathr2d2 (残宇)
看板Fortran
标题[问题] 能帮我看下程是哪里出问题吗??
时间Tue Apr 21 10:53:07 2009
写了一个用rk4普通解 解振动型态的二阶ode
可是结果却不如预期 结构周期不对
请问有没有高手能帮忙看看
结构方程式: 6000 * X'' + 622.08 * X' + 64499.72544 * X = - 6000 * XG
结构主频:1.0368
程式如下:
PROGRAM MAIN
IMPLICIT none
REAL*8 t, tstep, x(5)
REAL*8 XG
COMMON XG
INTEGER N, i, j, nsteps
OPEN(6, FILE='TIMEHIS.TXT')
N=2
nsteps=10000
tstep=0.002
x(1)=0.0
x(2)=0.0
DO 60 j = 1, nsteps
t = t+0.002
XG = SIN(t)
CALL rk41(t, x, tstep, N)
WRITE (6,*) t, X(1)
60 CONTINUE
CLOSE(6)
STOP
END
SUBROUTINE rk41(t, x, tstep, N)
IMPLICIT none
REAL*8 DERIV1, h, t, tstep, x(5)
REAL*8 k1(5), k2(5),k3(5), k4(5), temp1(5), temp2(5), temp3(5)
REAL*8 XG
COMMON XG
INTEGER i, N
h = tstep / 2
DO 10 i = 1,N
k1(i) = tstep * DERIV1(t, x, i)
temp1(i) = x(i) + 0.5*k1(i)
10 CONTINUE
DO 20 i = 1,N
k2(i) = tstep * DERIV1(t+h, temp1, i)
temp2(i) = x(i) + 0.5*k2(i)
20 CONTINUE
DO 30 i = 1,N
k3(i) = tstep * DERIV1(t+h, temp2, i)
temp3(i) = x(i) + k3(i)
30 CONTINUE
DO 40 i = 1,N
k4(i) = tstep * DERIV1(t + tstep, temp3, i)
x(i) = x(i) + (k1(i) + (2. * (k2(i) + k3(i))) + k4(i)) / 6.0
40 CONTINUE
RETURN
END
FUNCTION DERIV1(t, temp, i)
IMPLICIT none
REAL*8 DERIV1, t, temp(5)
REAL*8 XG
COMMON XG
INTEGER i
IF (i .EQ. 1) DERIV1 = temp(2)
IF (i .EQ. 2) DERIV1 = (-6000*XG - 622.08*temp(2) - 6449.72544*temp(1))/6000
RETURN
END
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.120.28.186
1F：→ YCTzeng:我测试过程式可以执行,所以问题应该不在程式。 04/26 23:35
2F：→ deathr2d2:嗯嗯 感谢您 我继续研究 04/27 15:12