NIAKRAIS
Newbie | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору Здравствуйте, я начинающий программист на фортране, написал программу интегрирования методом рунге-кутты 4 порядка системы уравнений DXj/Dt = (X(j+1) - X(j-2)) * X(j-1) - X(j) + 8.0. Проблема в том, что при большем времени(он же верний предел интегрирования) решение уходит на бесконечность, но известно что оно должно асцилировать возле 8.0. Есть подозрение, что бесконечность получается из-за погрешности в числах, хотя я использую REAL*8, чего в принципе должно хватать. Повышение точности до REAL*16 позволило несколько увеличить время, но проблему не решило. П.С. Тестировал на 1 уравнении, вес прекрасно работает. Код: IMPLICIT NONE INTEGER :: I,J REAL*8 :: H,T,TMAX, FUNC REAL*8, ALLOCATABLE :: X(:) J = 40 ALLOCATE( X(J) ) TMAX = 0.05D+00 H = 0.05D+00 T = 0.D+00 DO I = 1,J X(I) = 8.0D+00 END DO X(20) = X(20) * 1.001D+00 100 CONTINUE IF (T<TMAX) THEN CALL RUNGE(J,T,H,X) T=T+H WRITE(*,1000) (X(I),I=1,J) GOTO 100 END IF 1000 FORMAT(5F12.8) CLOSE(100) END ! J - КОЛИЧЕСТВО ТОЧЕК ! Т - ВРЕМЯ ! Х - КООРДИНАТА ! XCOF- ПРИБАВЛЯЮТСЯ К Х ПРИ ИНТЕГРИРОВАНИИ РУНГЕ-КУТТОЙ ! I - ТЕКУЩАЯ ТОЧКА REAL*8 FUNCTION FUNC(J,T,X,XCOF,I) INTEGER :: J,I REAL*8 :: T REAL*8 :: X(J) REAL*8 :: XCOF IF(I .EQ. J) FUNC = (X(1) - X(I-2)) * (X(I-1) + XCOF) - (X(I) + XCOF) + 8.D+0 IF(I-1 .EQ. 0) FUNC = (X(I+1) - X(J-1)) * (X(J) + XCOF) - (X(I) + XCOF) + 8.D+0 IF(I-2 .EQ. 0) FUNC = (X(I+1) - X(J)) * (X(I-1) + XCOF) - (X(I) + XCOF) + 8.D+0 IF((I>2) .AND. (I < J)) FUNC = (X(I+1) - X(I-2)) * (X(I-1) +XCOF) - (X(I) + XCOF) + 8.D+0 1000 FORMAT(5F12.8) END FUNCTION FUNC SUBROUTINE RUNGE(J,T,H,X) INTEGER, INTENT(IN) :: J REAL*8,INTENT(IN) :: T,H REAL*8,INTENT(INOUT) :: X(J) REAL*8 :: K1, K2, K3, K4, FUNC REAL*8, Allocatable :: TEMPX(:) INTEGER :: I ALLOCATE (TEMPX(J)) DO I=1,J K1 = FUNC(J, T, X, 0.0D+0,I) K2 = FUNC(J, T, X, (H/2.D+0) * K1,I) K3 = FUNC(J, T, X, (H/2.D+0) * K2,I) K4 = FUNC(J, T, X, H * K3,I) TEMPX(I) = X(I) + H/6.D+0 * (K1 + 2.D+0*K2 + 2.D+0*K3 + K4) END DO DO I=1,J X(I) = TEMPX(I) END DO DEALLOCATE (TEMPX) 1000 FORMAT(5F12.8) END SUBROUTINE RUNGE |
| Всего записей: 3 | Зарегистр. 24-11-2011 | Отправлено: 19:02 24-11-2011 | Исправлено: NIAKRAIS, 19:05 24-11-2011 |
|