[問題] 能幫我看下程是哪裡出問題嗎??

看板Fortran作者deathr2d2 (殘宇)時間15年前 (2009/04/21 10:53)推噓0(0推 0噓 2→)

留言2則, 2人參與討論串1/1

寫了一個用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