Константа в фортране

[Alexey_aka_Lave]

Код: [Выделить]


     program test
parameter n=10
real a
a=5
print *, a
endalexey@alexey-laptop:~$ gfortran test.for -o main1 -x fortran
test.for:3.72:

      real a                                                           
                                                                        1
ошибка: Unexpected data declaration statement at (1)

Код: [Выделить]


     program test
! parameter n=10
real a
a=5
print *, a
endalexey@alexey-laptop:~$ gfortran test.for -o main1 -x fortran
alexey@alexey-laptop:~$ ./main1
   5.0000000   
alexey@alexey-laptop:~$

как правильно описать константу, или что я не так делаю( На других компиляторах Фортрана все отлично)

[Kwah]

Код: [Выделить]


     program test
parameter n=10
real a
a=5
print *, a
end...
как правильно описать константу, или что я не так делаю( На других компиляторах Фортрана все отлично)

Взять n=10 в скобки.

Пользоваться implicit декларациями — моветон.

[Alexey_aka_Lave]

Огромнейшее спасибо))
Пользователь решил продолжить мысль 22 Мая 2009, 00:29:42:

Код: [Выделить]

REAL XP(DIMM), RP(DIMM), UP(DIMM)
REAL X0(DIMM)Вот такое описание не работает ругается.(в других компиляторах работает=))
Можно поподробней про:"Пользоваться implicit декларациями — моветон."

[Kwah]

Пожалуйста.

Код: [Выделить]

REAL XP(DIMM), RP(DIMM), UP(DIMM)
REAL X0(DIMM)Вот такое описание не работает ругается.(в других компиляторах работает=))

Что такое DIMM ? В общем, показывай целостные куски кода, в которых упоминаются по крайней мере декларации всего, что есть в проблемных строках и как именно ругается.

Другие компиляторы (интересно какие) думаю можно тоже настроить так, что они будут выдавать ошибки на твоём коде.

Можно поподробней про:"Пользоваться implicit декларациями — моветон."

Если ты в своей программе выше напишешь что-то типа

Код: [Выделить]

       program ...
       implicit none
       ...
то опять получишь ошибки. Прикол в том, что по умолчанию за именами переменных, начинающихся с конкретных букв, закреплён тип, implicit none это бредовое поведение по умолчанию отключает. Бредовое, ибо является ещё одинм источником ошибок. Приучайся писать, не опираясь на такие вот фишки.

[Alexey_aka_Lave]

(Нажмите, чтобы показать/скрыть)

   PROGRAM MAIN

   PARAMETER (DIMM=2)

   INTEGER DIMMENS
   REAL Kds, Tds, T01, T02, T3, T4
   REAL K01, K, F1, F2, FUNC1, FUNC2
   REAL EPS, RESULT(DIMM)
   REAL XP(DIMM), RP(DIMM), UP(DIMM)
   REAL X0(DIMM)
   COMMON /OBLX/ XP,RP,UP
   COMMON /K/ K
   COMMON /C/ Kds, Tds, T01, T02, K01

C   ïðéóáîéå ëïîóôáîô   
   
       DATA KDS,TDS,T01,T02,K01 /0.1,0.001,0.025,0.35,1.5/
       K=0.1

C   ***************
   
   EPS=0.1
   DIMMENS=DIMM

   T3=0.35
   T4=0.0085
C   F2=FUNC2(T3,T4)
   
   X0(1)=T3   
   X0(2)=T4
   TTT=F1(X0)
   PRINT *, TTT
   CALL OVRAZH(F1,EPS,DIMMENS,RESULT,X0)
   PRINT *, RESULT
   END

   REAL FUNCTION F1(X)   
   REAL X(2)
   F1=FUNC1(X(1),X(2))
   END

(Нажмите, чтобы показать/скрыть)

alexey@alexey-laptop:~$ gfortran main.for
main.for:8.23:

      REAL EPS, RESULT(DIMM)                                           
                       1
ошибка: Expression at (1) must be of INTEGER type
main.for:8.28:

      REAL EPS, RESULT(DIMM)                                           
                            1
ошибка: The module or main program array 'result' at (1) must have constant shape
main.for:9.24:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
                        1
ошибка: Expression at (1) must be of INTEGER type
main.for:9.29:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
                             1
ошибка: The module or main program array 'rp' at (1) must have constant shape
main.for:9.34:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
                                  1
ошибка: Expression at (1) must be of INTEGER type
main.for:9.39:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
                                       1
ошибка: The module or main program array 'up' at (1) must have constant shape
main.for:10.14:

      REAL X0(DIMM)                                                     
              1
ошибка: Expression at (1) must be of INTEGER type
main.for:10.19:

      REAL X0(DIMM)                                                     
                   1
ошибка: The module or main program array 'x0' at (1) must have constant shape
main.for:9.14:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
              1
ошибка: Expression at (1) must be of INTEGER type
main.for:9.19:

      REAL XP(DIMM), RP(DIMM), UP(DIMM)                                 
                   1
ошибка: The module or main program array 'xp' at (1) must have constant shape
alexey@alexey-laptop:~$

на OpenVMS (Не в курсе какой там компилятор))) отлично пашут)

[Kwah]

Мдя. Я же написал про IMPLICIT.

По умолчанию DIMM без явного объявления считается вещественным, ибо

IMPLICIT: default IMPLICIT REAL(A-H, O-Z), INTEGER(I-N)

(Нажмите, чтобы показать/скрыть)

Код: [Выделить]

      PROGRAM MAIN
      IMPLICIT NONE
      INTEGER DIMM

      PARAMETER (DIMM=2)
 
      INTEGER DIMMENS
      REAL Kds, Tds, T01, T02, T3, T4
      REAL K01, K, F1, F2, FUNC1, FUNC2
      REAL EPS, RESULT(DIMM)
      REAL XP(DIMM), RP(DIMM), UP(DIMM)
      REAL X0(DIMM)
      REAL TTT
      COMMON /OBLX/ XP,RP,UP
      COMMON /K/ K
      COMMON /C/ Kds, Tds, T01, T02, K01

C   ïðéóáîéå ëïîóôáîô   
   
       DATA KDS,TDS,T01,T02,K01 /0.1,0.001,0.025,0.35,1.5/
       K=0.1

C   ***************
   
      EPS=0.1
      DIMMENS=DIMM

      T3=0.35
      T4=0.0085
C   F2=FUNC2(T3,T4)
   
      X0(1)=T3   
      X0(2)=T4
      TTT=F1(X0)
      PRINT *, TTT
      CALL OVRAZH(F1,EPS,DIMMENS,RESULT,X0)
      PRINT *, RESULT
      END

      REAL FUNCTION F1(X)   
      REAL X(2)
      F1=FUNC1(X(1),X(2))
      END



[Alexey_aka_Lave]

Знал я про INTEGER(I-N) тока сразу не понял что про них родимых речь идет))
теперь вообще супер))
ЗЫ а как константу объявить определенного типа, пока добавил IDIMM, так на будущее что бы знать ?(раньше не парился с этим))

[Kwah]

Заглянуть под спойлер. Там написано.

[Alexey_aka_Lave]

не заметил сначала, все супер отлично пашет.
Вот тока осталось с дебагером разобраться). GDB конечно хорошо, но мне бы попроще и по нагляднее.
Хотелось что ни-будь как на OpenVMS
Пол экрана код, указатель показывает выполняемую строку, можно проверить значение всех переменных

[Kwah]

К gdb есть несколько граф.интерфейсов.

[Alexey_aka_Lave]

ddd весьма неплох тока ещё не полностью разобрался.

(Нажмите, чтобы показать/скрыть)

Код: [Выделить]

PROGRAM pp
IMPLICIT NONE

INTEGER iwork(5),n,iflag,i

parameter (n=4)

real x(n),work(6*n+3),u,t1,t2,t3,t4,tds,k01,k02,kds,xx,model,sum

real xt,x1st,xst,to,tk,y,t,re,s,k,flag,ae,est,e,dt

        external fun

common/u/ u,t1,t2,t3,t4,tds,k01,k02,kds

flag=0        

t1=0.025

        t2=0.35

        t3=0.35

        t4=0.0085

        tds=0.001





        k01=1.5

        k02=0.7

        kds=0.1  

dt=0.0001

        T=0


xst=0
x1st=0


k=1

to=t

DO 11 I=1,N

11 x(i)=0.0

        tk=6.28

re=1.e-6

ae=0.



iflag=1

30      continue

75 CALL RKF45(fun,n,x,t,to,re,ae,iflag,work,iwork)

go to(88,82,83,84,85,86,87,88),iflag

83 iflag=2

re=re*10

go to 75

84      iflag=2

go to 75

85 iflag=2

ae=1.e-6

go to 75

86 iflag=2

re=re*10

go to 75

87 iflag=2

go to 75

88 stop111

82      T=TO

y=1.05*sin(0.5*t)
est=xst-kds*x1st

        e=y-Kds*x(4)

        S=S+e*e*dt+(k*(e-est)*(e-est)/dt)

TO=TO+DT

        RE=1.E-6

        AE=0.


! WRITE (6,*)T          !!!!!!!!!!!!!!!!!!!!!!!!!!!!
XST=y
X1ST=X(4)

est=e
        IF(T.LT.TK)GOTO 30
WRITE (2,*)T
print *,s
close(2)

end



        SUBROUTINE FUN(t,a,ap)

        REAL a(4),ap(4),xx,t1,t2,tds,t3,t4,k01,k02,kds,u,w,alp,c,b,y
        common/u/ u,t1,t2,t3,t4,tds,k01,k02,kds

w=2*sin(50*t)

xx=1.05*sin(0.5*t)

b=xx-a(1)

alp=(t3/t4)*b

ap(2)=b/t4

if ((a(2).GE.10).and.(ap(2).GT.0)) ap(2)=0

if ((a(2).LE.-10).and.(ap(2).GT.0)) ap(2)=0

ap(1)=(a(4)*kds-a(1))/tds

u=a(2)+alp

u=(abs(u+10)-abs(u-10))/2

ap(3)=(k01*u-a(3))/t1

c=a(3)-w

ap(4)=(c*k02-a(4))/t2




RETURN

        END





Если раскоментировать !WRITE (6,*)T программа зависает отказывается работать .
Без печати работает вроде нормально. В чем проблема кто подскажет.

[Kwah]

Файл-то, которому присвоен нумер 6 открыт? На запись?

[Alexey_aka_Lave]

Это вывод на экран.
если написать в этом месте print *,T  происходит тоже самое.

[Kwah]

Это вывод на экран.
если написать в этом месте print *,T  происходит тоже самое.

Ну, тебе видней, потому что из приведённого кода этого понять невозможно. В общем, показывай код, или по крайней мере рассказывай где и каким образом открываются файлы на запись... А вообще, кроме этого в данной строке не видно ничего криминального, ошибка скорее всего в другом месте, а здесь она лишь проявляет себя.

ЗЫ. За подобный стиль кодирования нужно руки отрывать.

[Alexey_aka_Lave]

Этот же код на ОпенВмс на экран (и в файл) отлично печатает

(Нажмите, чтобы показать/скрыть)

Код: [Выделить]

 SUBROUTINE RKF45(F,NEQN,Y,T,TOUT,RELERR,ABSERR,
     *                 IFLAG,WORK,IWORK)
C
C     METOД PУHГE-KУTTA ФEЛЬБEPГA ЧETBEPTOГO-ПЯTOГO ПOPЯДKA
C
C     COCTABИTEЛИ ПPOГPAMMЫ-H.A.WATTS,L.F.SHAMPINE
C                SANDIA LABORATORIES
C              ALBUQUERQUE, NEW MEXICO
C
C     RKF45 ПPEДHAЗHAЧEHA ГЛABHЫM OБPAЗOM ДЛЯ PEШEHИЯ
C     HEЖECTKИX И CЛAБO ЖECTKИX ДИФФEPEHЦИAЛЬHЫX УPABHEHИЙ,
C     KOГДA BЫЧИCЛEHИE ПPOИЗBOДHЫX HE CЛИШKOM ДOPOГOCTOЯЩEE.
C     RKF45, BOOБЩE ГOBOPЯ,HE CЛEДУET ИCПOЛЬЗOBATЬ
C     ECЛИ ПOЛЬЗOBATEЛЮ TPEБУETCЯ BЫCOKAЯ TOЧHOCTЬ
C
C     PEЗЮME
C
C     ПOДПPOГPAMMA RKF45 ИHTEГPИPУET CИCTEMУ ИЗ NEQN
C     OБЫKHOBEHHЫX ДИФФEPEHЦИAЛЬHЫX УPABHEHИЙ ПEPBOГO
C     ПOPЯДKA CЛEДУЮЩEГO BИДA:
C
C            DY(I)/DT=F(T,Y(1),Y(2),...,Y(NEQN),
C
C     ГДE Y(I) ЗAДAHЫ B T.
C     OБЫЧHO ПOДПPOГPAMMУ ПPИMEHЯЮT ДЛЯ ИHTEГPИPOBAHИЯ
C     OT T ДO TOUT, OДHAKO EE MOЖHO ИCПOЛЬЗOBATЬ И KAK
C     OДHOШAГOBЫЙ ИHTEГPATOP,ЧTOБЫ ПPOДOЛЖИTЬ PEШEHИE HA
C     OДИH ШAГ B HAПPABЛEHИИ TOUT.HA BЫXOДE ПAPAMETPAM,
C     ФИГУPИPУЮЩИM B CПИCKE BЫЗOBA, ПPИCBAИBAЮTCЯ ЗHAЧEHИЯ,
C     HEOБXOДИMЫE ДЛЯ ПPOДOЛЖEHИЯ ИHTEГPИPOBAHИЯ. ПOЛЬЗO-
C     BATEЛЮ HУЖHO ЛИШЬ EЩE PAЗ OБPATИTЬCЯ K RKF45
C     (И,BOЗMOЖHO,OПPEДEЛИTЬ HOBOE ЗHAЧEHИE ДЛЯ TOUT).
C     B ДEЙCTBИTEЛЬHOCTИ RKF45-ЭTO ПPOГPAMMA ИHTEPФEЙCA,
C     KOTOPAЯ BЫЗЫBAET ПOДПPOГPAMMУ RKFS, OCУЩECTBЛЯЮЩУЮ
C     ПPOЦECC PEШEHИЯ.RKFS B CBOЮ OЧEPEДЬ BЫЗЫBAET ПOДПPOГ-
C     PAMMУ FEHL, KOTOPAЯ BЫЧИCЛЯET ПPИБЛИЖEHHOE PEШEHИE
C     HA OДИH ШAГ.
C
C     RKF45 ИCПOЛЬЗУET METOД PУHГE-KУTTA-ФEЛЬБEPГA, OПИCAHHЫЙ
C     B CЛEДУЮЩEЙ ПУБЛИKAЦИИ:E.FEHLBERG,LOW-ORDER CLASSICAL
C     RUNGE-KUTTA FORMULAS WITH STEPSIZE CONTROL,NASA TR R-315
C
C     CTИЛЬ PAБOTЫ ПPOГPAMMЫ RKF45 ИЛЛЮCTPИPУETCЯ B CЛEДУЮЩИX
C     ПУБЛИKAЦИЯX:L.F.SHAMPINE,H.A.WATTS,S.DAVENPORT, SOLVING
C     NON-STIFF ORDINARY DIFFERENTIAL EQUATIONS-THE STAT OF
C     THE ART,SANDIA LABORATORIES REPORT SAND75-0182,SIAM
C     REVIEW,18(1976), N3,376-411.
C
C     ПAPAMETPЫ ПPOГPAMMЫ:
C
C     F       -ПOДПPOГPAMMA F(T,Y,YP) ДЛЯ BЫЧИCЛEHИЯ
C              ПPOИЗBOДHЫX YP(I)=DY(I)/DT
C     NEQN    -ЧИCЛO ИHTEГPИPУEMЫX УPABHEHИЙ
C     Y(*)    -PEШEHИE B TOЧKE T
C     T       -HEЗABИCИMAЯ ПEPEMEHHAЯ
C     TOUT    -TOЧKA BЫXOДA,B KOTOPOЙ HУЖHO OПPEДEЛИTЬ
C              ЗHAЧEHИE PEШEHИЯ
C     RELERR  -ГPAHИЦA OTHOCИTEЛЬHOЙ ПOГPEШHOCTИ
C              ДЛЯ TECTA ЛOKAЛЬHOЙ OШИБKИ.
C     ABSERR  -ГPAHИЦA ABCOЛЮTHOЙ ПOГPEШHOCTИ
C              ДЛЯ TECTA ЛOKAЛЬHOЙ OШИБKИ.
C              HA KAЖДOM ШAГE ПPOГPAMMA TPEБУET BЫПOЛHEHИЯ УCЛOBИЯ
C              ABS(LOCAL ERROR).LE.RELERR*ABS(Y)+ABSERR
C              ДЛЯ KAЖДOЙ KOMПOHEHTЫ BEKTOPOB ЛOKAЛЬHOЙ
C              OШИБKИ И PEШEHИЯ
C     IFLAG   -УKAЗATEЛЬ PEЖИMA ИHTEГPИPOBAHИЯ.
C     WORK(*) -MACCИB,COДEPЖAЩИЙ ИHФOPMAЦИЮ,BHУTPEHHЮЮ ДЛЯ RKF45,
C              KOTOPAЯ HEOБXOДИMA ПPИ ПOCЛEДУЮЩИX BЫЗOBAX.EГO
C              PAЗMEPHOCTЬ ДOЛЖHA БЫTЬ HE MEHЬШE 3+6*NEQN
C     IWORK(*)-ЦEЛЫЙ MACCИB,COДEPЖAЩИЙ ИHФOPMAЦИЮ,BHУTPEHHЮЮ ДЛЯ
C              RKF45,KOTOPAЯ HEOБXOДИMA ПPИ ПOCЛEДУЮЩИX BЫЗOBAX.
C              EГO PAЗMEPHOCTЬ ДOЛЖHA БЫTЬ HE MEHЬШE 5.
C
C     ПEPBOE OБPAЩEHИE K RKF45
C
C     ПOЛЬЗOBATEЛЬ ДOЛЖEH ПPEДУCMOTPETЬ B CBOEЙ BЫЗЫBAЮЩEЙ
C     ПPOГPAMME ПAMЯTЬ ДЛЯ CЛEДУЮЩИX MACCИBOB, ФИГУPИPУЮЩИX
C     B CПИCKE BЫЗOBA- Y(NEQN), WORK(3+6*NEQN), IWORK(5);
C     KPOME TOГO, OH ДOЛЖEH OБ'ЯBИTЬ F B OПEPATOPE EXTERNAL,
C     ПOДГOTOBИTЬ ПOДПPOГPAMMУ F(T,Y,YP) И ПPИCBOИTЬ HAЧAЛЬ-
C     HЫE ЗHAЧEHИЯ ПAPAMETPAM-
C
C     NEQN  -ЧИCЛO ИHTEГPИPУEMЫX УPABHEHИЙ (NEQN.GE.1)
C     Y(*)  -BEKTOP HAЧAЛЬHЫX УCЛOBИЙ
C     T     -HAЧAЛЬHAЯ TOЧKA ИHTEГPИPOBAHИЯ,
C            T ДOЛЖHO БЫTЬ ПEPEMEHHOЙ.
C     TOUT  -TOЧKA BЫXOДA,B KOTOPOЙ HУЖHO HAЙTИ ЗHAЧEHИE
C            PEШEHИЯ. T=TOUT BOЗMOЖHO ЛИШЬ ПPИ ПEPBOM
C            OБPAЩEHИИ.B ЭTOM CЛУЧAE BЫXOД ИЗ RKF45 ПPOИ-
C            CXOДИT CO ЗHAЧEHИEM ПAPAMETPA IFLAG=2,ECЛИ
C            MOЖHO ПPOДOЛЖATЬ ИHTEГPИPOBAHИE.
C     RELERR-ГPAHИЦA ДЛЯ OTHOCИTEЛЬHOЙ ЛOKAЛЬHЫЙ ПOГPEШHOCTEИ.
C     ABSERR-ГPAHИЦA ДЛЯ AБCOЛЮTHOЙ    ЛOKAЛЬHЫЙ ПOГPEШHOCTEИ.
C            ЭTИ ГPAHИЦЫ ДOЛЖHЫ БЫTЬ HEOTPИЦATEЛЬHЫ.
C            RELERR ДOЛЖHA БЫTЬ ПEPEMEHHOЙ,A ABSERR MOЖET
C            БЫTЬ И KOHCTAHTOЙ.ПPOГPAMME, BOOБЩE ГOBOPЯ
C            HE CЛEДУET ЗAДABATЬ ГPAHИЦУ ДЛЯ OTHOCИTEЛЬHOЙ
C            OШИБKИ,MEHЬШУЮ, ЧEM ПPИMEPHO 1.E-7. ДAБЫ ИЗБEЖATЬ
C            TPУДHOCTEЙ ,CBЯЗAHHЫX C OЧEHЬ BЫCOKИMИ ЗAПPOCAMИ
C            K TOЧHOCTИ, ПPOГPAMMA TPEБУET,ЧTOБЫ RELERR
C            БЫЛA БOЛЬШE, ЧEM HEKOTOPЫЙ ПAPAMETP OTHOCИTEЛЬHOЙ
C            OШИБKИ,BЫЧИCЛЯEMЫЙ BHУTPИ EE И ЗABИCЯЩИЙ OT
C            MAШИHЫ.B ЧACTHOCTИ,HE PAЗPEШAETCЯ ЗAДAHИE TOЛЬKO
C            AБCOЛЮTHOЙ OШИБKИ.ECЛИ ЖE ЗAДAHO ЗHAЧEHИE RELERR,
C            MEHЬШEE ДOПУCTИMOГO, TO RKF45 УBEЛИЧИBAET RELERR
C            HAДЛEЖAЩИM OБPAЗOM И BOЗBPAЩAET УПPABЛEHИE ПOЛЬ-
C            ЗOBATEЛЮ, ПPEЖДE ЧEM ПPOДOЛЖATЬ ИHTEГPИPOBAHИE.
C     IFLAG-+1,-1.ЭTO УKAЗATEЛЬ HACTPOЙKИ ПPOГPAMMЫ ДЛЯ KAЖДOЙ
C            HOBOЙ ЗAДAЧИ. HOPMAЛЬHOE BXOДHOE ЗHAЧEHИE PABHO+1.
C            ПOЛЬЗOBATEЛЬ ДOЛЖEH ЗAДABATЬ IFLAG=-1 ЛИШЬ B TOM
C            CЛУЧAE,  KOГДA HEOБXOДИMO УПPABЛEHИE OДHOШAГOBЫM
C            ИHTEГPATOPOM.B ЭTOM CЛУЧAE RKF45 ПЫTAETCЯ ПPOДOЛЖИTЬ
C            PEШEHИE HA OДИH ШAГ B HAПPABЛEHИИ TOUT ПPИ KAЖДOM
C            OЧEPEДHOM BЫЗOBE. ПOCKOЛЬKУ ЭTOT PEЖИM PAБOTЫ
C            BECЬMA HEЭKOHOMИЧEH, EГO CЛEДУET ПPИMEHЯTЬ
C            ЛИШЬ B CЛУЧAE KPAЙHEЙ HEOБXOДИMOCTИ.
C
C     ИHФOPMAЦИЯ HA BЫXOДE
C
C     Y(*)    -PEШEHИE B TOЧKE T
C     T       -ПOCЛEДHЯЯ TOЧKA,ДOCTИГHУTAЯ ПPИ ИHTEГPИPOBAHИИ.
C     IFLAG=2 -ПPИИHTEГPИPOBAHИИ ДOCTИГHУTO TOUT.ЭTO ЗHAЧEHИE
C              ПAPAMETPA УKAЗЫBAET HA УCПEШHЫЙ BЫXOД И
C              ЯBЛЯETCЯ HOPMAЛЬHЫM PEЖИMOM ДЛЯ ПPOДOЛЖEHИЯ
C              ИHTEГPИPOBAHИЯ.
C          =3 -ИHTEГPИPOBAHИE HE БЫЛO ЗAKOHЧEHO ИЗ-ЗA TOГO,
C              ЧTO ЗAДAHHOE ЗHAЧEHИE ГPAHИЦЫ ДЛЯ OTHOCИTEЛЬHOЙ
C              OШИБKИ OKAЗAЛOCЬ CЛИШKOM MAЛO. ДЛЯ ПPOДOЛЖEHИЯ
C              ИHTEГPИPOBAHИЯ RELERR БЫЛO HAДЛEЖAЩИM OБPAЗOM
C              УBEЛИЧEHO.
C          =4 -ИHTEГPИPOBAHИE HE БЫЛO ЗAKOHЧEHO ИЗ-ЗA TOГO,
C              ЧTO ПOTPEБOBAЛOCЬ БOЛEE 3000 BЫЧИCЛEHИЙ ПPO-
C              ИЗBOДHOЙ.ЭTO COOTBETCTBYET ПPИБЛИЗИTEЛЬHO
C              500 ШAГAM.
C          =5 -ИHTEГPИPOBAHИE HE БЫЛO ЗAKOHЧEHO ИЗ-ЗA TOГO,
C              ЧTO PEШEHИE OБPATИЛOCЬ B HYЛЬ,BCЛEДCTBИE ЧEГO
C              TECT TOЛЬKO OTHOCИTEЛЬHOЙ OШИБKИ HE ПPOXOДИT.
C              ДЛЯ ПPOДOЛЖEHИЯ HEOБXOДИMO HEHYЛEBOE ЗHAЧEHИE
C              ПAPAMETPA ABSERR. ИCПOЛЬЗOBAHИE HA OДИH ШAГ
C              PEЖИMA ПOШAГOBOГO ИHTEГPИPOBAHИЯ ЯBЛЯETCЯ
C              PAЗYMHЫM BЫXOДOM ИЗ ПOЛOЖEHИЯ.
C          =6 -ИHTEГPИPOBAHИE  HE БЫЛO ЗAKOHЧEHO ИЗ-ЗA TOГO,
C              ЧTO TPEБYEMAЯ TOЧHOCTЬ HE MOГЛA БЫTЬ ДOCTИГHУTA
C              ДAЖE ПPИ HAИMEHЬШEЙ ДOПУCTИ MOЙ BEЛИЧИHE ШAГA.
C              ПOЛЬЗOBATEЛЬ ДOЛЖEH УBEЛИЧИTЬ ГPAHИЦУ ПOГPEШ-
C              HOCTИ,ПPEЖДE ЧEM MOЖHO БУДET ПOПЫTATЬCЯ
C              ПPOДOЛЖATЬ ИHTEГPИPOBAHИE.
C          =7 -ПO BCEЙ BИДИMOCTИ, RKF45 HEЭФФEKTИBHA ПPИ
C              PEШEHИИ ЭTOЙ ЗAДAЧИ. CЛИШKOM БOЛЬШOE ЧИCЛO
C              TPEБУEMЫX BЫXOДHЫX TOЧEK ПPEПЯTCTBУET BЫБOPУ
C              ECTECTBEHHOЙ BEЛИЧИHЫ ШAГA.CЛEДУET ИCПOЛЬЗOBATЬ
C              PEЖИM ПOШAГOBOГO ИHTEГPИPOBAHИЯ.
C          =8 -HEПPABИЛЬHOE ЗAДAHИE BXOДHЫX ПAPAMETPOB.ЭTO
C              ЗHAЧEHИE ПOЯBЛЯETCЯ,ECЛИ ДOПУЩEHA OДHA ИЗ
C              CЛEДУЮЩИX OШИБOK-
C                              NEQN.LE.0
C                  T=TOUT  И  IFLAG.NE.+1  ИЛИ -1
C                  RELERR  ИЛИ  ABSERR.LT.0
C                  IFLAG.EQ.0  ИЛИ .LT.-2  ИЛИ .GT.8
C     WORK(*) -ИHФOPMAЦИЯ, KOTOPAЯ OБЫЧHO HE ПPEДCTABЛЯET ИHTE-
C              PECA ДЛЯ ПOЛЬЗOBATEЛЯ, HO HEOБXOДИMA ПPИ ПOCЛE-
C              ДУЮЩИX BЫЗOBAX. WORK(1),...,WORK(NEQN) COДEPЖAT
C              ПEPBЫE ПPOИЗBOДHЫE BEKTOPA PEШEHИЯ Y B TOЧKE T.
C              WORK(NEQN+1) XPAHИT BEЛИЧИHУ ШAГA H,C KOTOPOЙ
C              MOЖHO ПOПЫTATЬCЯ ПPOBECTИ CЛEДУЮЩИЙ ШAГ.
C     IWORK(*) -ИHФOPMAЦИЯ, KOTOPAЯ OБЫЧHO HE ПPEДCTABЛЯET ИHTE-
C               PECA ДЛЯ ПOЛЬЗOBATEЛЯ, HO HEOБXOДИMA ПPИ ПOCЛE-
C               ДУЮЩИX BЫЗOBAX. B IWORK(1) COДEPЖИTCЯ
C               CЧETЧИK ЧИCЛA BЫЧИCЛEHИЙ ПPOИЗBOДHЫX.
C
C     ПOCЛEДУЮЩИE OБPAЩEHИЯ K RKF45
C
C          HA BЫXOДE ПOДПPOГPAMMЫ RKF45 ИMEETCЯ BCЯ ИHФOPMAЦИЯ,
C     HEOБXOДИMAЯ  ДЛЯ ПPOДOЛЖEHИЯ ИHTEГPИPOBAHИЯ.ECЛИ ПPИ
C     ИHTEГPИPOBAHИИ ДOCTИГHУTO TOUT,TO ПOЛЬЗOBATEЛЮ ДOCTA-
C     TOЧHO OПPEДEЛИTЬ HOBOE ЗHAЧEHИE  TOUT И CHOBA OБPATИTЬ-
C     CЯ K RKF45.
C          B PEЖИME ПOШAГOBOГO ИHTEГPИPOBAHИЯ (IFLAG=-2)
C     ПOЛЬЗOBATEЛЬ ДOЛЖEH ИMETЬ B BИДУ,ЧTO KAЖДЫЙ ШAГ
C     BЫПOЛHЯETCЯ B HAПPABЛEHИИ TEKУЩEГO ЗHAЧEHИЯ TOUT
C     (CИГHAЛИЗИPУEMOM  ИЗMEHEHИEM IFLAG HA 2). ПOЛЬЗOBATEЛЬ
C     ДOЛЖEH ЗAДATЬ HOBOE ЗHAЧEHИE TOUT И ПEPEOПPEДEЛИTЬ
C     IFLAG HA -2, ЧTOБЫ ПPOДOЛЖATЬ B PEЖИME ПOШAГOBOГO
C     ИHTEГPИPOBAHИЯ.
C          ECЛИ ИHTEГPИPOBAHИE HE БЫЛO ЗAKOHЧEHO,HO
C     ПOЛЬЗOBATEЛЬ XOЧET ПPOДOЛЖATЬ (CЛУЧAИ IFLAG=3,4), OH
C     ПOПPOCTУ CHOBA OБPAЩAETCЯ K RKF45.ПPИ IFLAG=3 ПAPA-
C     METP RELERR БЫЛ  ИЗMEHEH HAДЛEЖAЩИM ДЛЯ ПPOДOЛЖEHИЯ
C     ИHTEГPИPOBAHИЯ OБPAЗOM.B CЛУЧAE IFLAG=4 CЧETЧИK
C     ЧИCЛA ЗHAЧEHИЙ ФУHKЦИИ БУДET ПEPEOПPEДEЛEH HA 0, И
C     БУДУT PAЗPEШEHЫ EЩE 3000 BЫЧИCЛEHИЙ ФУHKЦИИ.
C          OДHAKO B CЛУЧAE IFLAG=5, ПPEЖДE ЧEM MOЖHO БУДET
C     ПPOДOЛЖATЬ ИHTEГPИPOBAHИE,ПOЛЬЗOBATEЛЬ ДOЛЖEH CHAЧAЛA
C     ИЗMEHИTЬ KPИTEPИЙ OШИБKИ, ЗAДAB ПOЛOЖИTEЛЬHOE ЗHAЧEHИE
C     ДЛЯ ABSERR. ECЛИ OH HE CДEЛAET ЭTOГO, BЫПOЛHEHИE ПPO-
C     ГPAMMЫ БУДET ПPEKPAЩEHO.
C        TOЧHO TAK ЖE,B CЛУЧAE IFLAG=6,ПPEЖДE ЧEM ПPOДOЛ-
C     ЖATЬ ИHTEГPИPOBAHИE,ПOЛЬЗOBATEЛЮ HEOБXOДИMO ПEPEOПPE-
C     ДEЛИTЬ IFLAG HA 2 (ИЛИ -2, ECЛИ ИCПOЛЬЗУETCЯ PEЖИM
C     ПOШAГOBOГO ИHTEГPИPOBAHИЯ) И УBEЛИЧИTЬ ЗHAЧEHИE ДЛЯ
C     ABSERR ЛИБO RELERR,ЛИБO И ДЛЯ TOГO,И ДЛЯ ДPУГOГO.
C     ECЛИ ЭTO HE БУДET CДEЛAHO,BЫПOЛHEHИE ПPOГPAMMЫ
C     ПPEKPAЩAETCЯ. ПOЯBЛEHИE IFLAG=6 УKAЗЫBAET HA HEPEГУ-
C     ЛЯPHOCTЬ (PEШEHИE БЫCTPO MEHЯETCЯ ИЛИ, BOЗMOЖHO,
C     ИMEETCЯ OCOБEHHOCTЬ),И ЧACTO B ПOДOБHЫX CЛУЧAЯX
C     HE ИMEET CMЫCЛA ПPOДOЛЖATЬ ИHTEГPИPOBAHИE.
C          ECЛИ БУДET ПOЛУЧEHO ЗHAЧEHИE IFLAG=7,TO ПOЛЬЗO-
C     BATEЛЬ ДOЛЖEH ПEPEЙTИ K PEЖИMУ ПOШAГOBOГO ИHTEГPИPO-
C     BAHИЯ C BEЛИЧИHOЙ ШAГA,OПPEДEЛЯEMOЙ ПPOГPAMMOЙ, ИЛИ
C     PACCMOTPETЬ BOЗMOЖHOCTTЬ ПEPEXOДA HA ПPOГPAMMЫ METOДOB
C     AДAMCA.ECЛИ BCE ЖE ПOЛЬЗOBATEЛЬ XOЧET ПPOДOЛЖATЬ
C     ИHTEГPИPOBAHИE ПO ПOДПPOГPAMME RKF45,OH ДOЛЖEH ДO HOBOГO
C     OБPAЩEHИЯ K HEЙ ПEPEOПPEДEЛИTЬ IFLAG HA 2.B ПPOTИBHOM
C     CЛУЧAE BЫПOЛHEHИE ПPOГPAMMЫ БУДET ПPEKPAЩEHO.
C          ECЛИ ПOЛУЧEHO ЗHAЧEHИE IFLAG=8,TO ИHTEГPИPOBAHИE
C     HEЛЬЗЯ ПPOДOЛЖATЬ,ПOKA HE БУДУT ИCПPABЛEHЫ OШИБOЧHЫE
C     BXOДHЫE ПAPAMETPЫ. HУЖHO OTMETИTЬ,ЧTO MACCИBЫ WORK И
C     IWORK COДEPЖAT ИHФOPMAЦИЮ,HEOБXOДИMУЮ ДЛЯ ДAЛЬHEЙШEГO
C     ИHTEГPИPOBAHИЯ.ПOЭTOMУ B ЭTИ MACCИBЫ HEЛЬЗЯ BHOCИTЬ
C     ИЗMEHEHИЙ.
C
      EXTERNAL F
      INTEGER NEQN,IFLAG,IWORK(5)
      REAL Y(NEQN),T,TOUT,RELERR,ABSERR,WORK(1)
C
C     ECЛИ TPAHCЛЯTOP ПPOBEPЯET ИHДEKCЫ, TO ЗAMEHИTЬ
C     WORK(1) HA WORK(3+6*NEQN)
C
      INTEGER K1,K2,K3,K4,K5,K6,K1M
C
C     BЫЧИCЛИTЬ ИHДEKCЫ ДЛЯ PACЩEПЛEHИЯ PAБOЧEГO MACCИBA
C
      K1M=NEQN+1
      K1=K1M+1
      K2=K1+NEQN
      K3=K2+NEQN
      K4=K3+NEQN
      K5=K4+NEQN
      K6=K5+NEQN
C
C     ЭTA ПPOMEЖYTOЧHAЯ ПPOГPAMMA ПPOCTO COKPAЩAET ДЛЯ
C     ПOЛЬЗOBATEЛЯ ДЛИHHЫЙ CПИCOK BЫЗOBA ПYTEM PACЩEПЛEHИЯ
C     ДBYX PAБOЧИX MACCИBOB. ECЛИ ЭTO HE COBMECTИMO C
C     TPAHCЛЯTOPOM,TO OH ДOЛЖEH OБPAЩATЬCЯ HEПOCPEДCTBEHHO
C     K ПOДПPOГPAMME RKFS .
C
      CALL RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,
     *          WORK(1),WORK(K1M),WORK(K1),WORK(K2),
     *          WORK(K3),WORK(K4),WORK(K5),WORK(K6),
     *          WORK(K6+1),IWORK(1),IWORK(2),
     *          IWORK(3),IWORK(4),IWORK(5))
      RETURN
      END




(Нажмите, чтобы показать/скрыть)

Код: [Выделить]

 SUBROUTINE RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,
     *                YP,H,F1,F2,F3,F4,F5,SAVRE,SAVAE,
     *                NFE,KOP,INIT,JFLAG,KFLAG)
C
C     METOД PУHГE-KУTTA-ФEЛЬБEPГA ЧETBEPTOГO-ПЯTOГO ПOPЯДKA
C     RKFS ИHTEГPИPУET CИCTEMУ OБЫKHOBEHHЫX ДИФФE-
C     PEHЦИAЛЬHЫX УPABHEHИЙ ПEPBOГO ПOPЯДKA(CM. KOM-
C     MEHTAPИЙ K RKF45). MACCИBЫ YP,F1,F2,F3,F4 И F5
C     (PAЗMEPHOCTИ ПO KPAЙHEЙ MEPE NEQN) И ПEPEMEH-
C     HЫE H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG И KFLAG
C     ИCПOЛЬЗУЮTCЯ BHУTPИ ПPOГPAMMЫ И BЫHECEHЫ B
C     CПИCOK BЫЗOBA,ЧTOБЫ COXPAHИTЬ ИX OПPEДEЛEH-
C     HOCTЬ ПPИ ПOBTOPHOM OБPAЩEHИИ.ПOЭTOMУ ИX ЗHA-
C     ЧEHИЯ HE ДOЛЖHЫ  ИЗMEHЯTЬCЯ  ПOЛЬЗOBATEЛEM.
C     BOЗMOЖHЫЙ ИHTEPEC ПPEДCTABЛЯЮT ПAPAMETPЫ
C     YP  -ПPOИЗBOДHAЯ BEKTOPA PEШEHИЯ B TOЧKE T
C     H   -ПPEДПOЛAГAEMЫЙ PAЗMEP ШAГA ДЛЯ OЧEPEДHOГO ЭTAПA
C     NFE -CЧETЧИK ЧИCЛA BЫЧИCЛEHИЙ ФУHKЦИИ
C
      LOGICAL HFAILD,OUTPUT
C
      INTEGER NEQN,IFLAG,NFE,KOP,INIT,JFLAG,KFLAG
      REAL Y(NEQN),T,TOUT,RELERR,ABSERR,H,YP(NEQN),
     *     F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN),
     *     SAVRE,SAVAE
C
C     EXTERNAL F
C
      REAL A,AE,DT,EE,EEOET,ESTTOL,ET,HMIN,REMIN,
     *     RER,S,SCAL,TOL,TOLN,U26,EPSP1,EPS,YPK
C
      INTEGER K,MAXNFE,MFLAG
C
      REAL AMAX1,AMIN1
C
C     REMIN-ЭTO MИHИMAЛЬHOE ДOПУCTИMOE ЗHAЧEHИE ДЛЯ
C     RELERR.ПOПЫTKИ ПOЛУЧИTЬ ПO ЭTOЙ ПOДПPOГPAMME
C     БOЛEE BЫCOKУЮ TOЧHOCTЬ OБЫЧHO CTOЯT OЧEHЬ
C     ДOPOГO И ЗAЧACTУЮ БEЗУCПEШHЫ.
C
      DATA REMIN/1.E-12/
C
C     CTOИMOCTЬ CЧETA KOHTPOЛИPУETCЯ TPEБOBAHИEM,
C     ЧTOБЫ KOЛИЧECTBO BЫЧИCЛEHИЙ ФУHKЦИИ БЫЛO OГ-
C     PAHИЧEHO BEЛИЧИHOЙ,ПPИБЛИЗИTEЛЬHO PABHOЙ ЗHA-
C     ЧEHИЮ ПAPAMETPA MAXNFE.ПPИHЯTOE ЗДECЬ ЗHAЧE-
C     HИE ПPИMEPHO COOTBETCTBУET 500 ШAГAM.
C
      DATA MAXNFE/3000/
C
C     ПPOBEPИTЬ BXOДHЫE ПAPAMETPЫ
C
      IF(NEQN.LT.1)GO TO 10
      IF((RELERR.LT.0.0).OR.(ABSERR.LT.0.0))GO TO 10
      MFLAG=IABS(IFLAG)
      IF((MFLAG.EQ.0).OR.(MFLAG.GT.8))GO TO 10
      IF(MFLAG.NE.1)GO TO 20
C
C     ПEPBЫЙ BЫЗOB,BЫЧИCЛИTЬ MAШИHHOE ЭПCИЛOH
C
      EPS=1.0
    5 EPS=EPS/2.0
      EPSP1=EPS+1.
      IF(EPSP1.GT.1.)GO TO 5
      U26=26.*EPS
      GO TO 50
C
C     OШИБKA BXOДHOЙ ИHФOPMAЦИИ
C
   10 IFLAG=8
      RETURN
C
C     ПPOBEPИTЬ BOЗMOЖHOCTИ ПPOДOЛЖEHИЯ
C
   20 IF((T.EQ.TOUT).AND.(KFLAG.NE.3))GO TO 10
      IF(MFLAG.NE.2)GO TO 25
C
C     IFLAG=+2 ИЛИ -2
C
      IF((KFLAG.EQ.3).OR.(INIT.EQ.0))GO TO 45
      IF(KFLAG.EQ.4)GO TO 40
      IF((KFLAG.EQ.5).AND.(ABSERR.EQ.0.0))GO TO 30
      IF((KFLAG.EQ.6).AND.(RELERR.LE.SAVRE).AND.
     *(ABSERR.LE.SAVAE))GO TO 30
      GO TO 50
C
C     IFLAG=3,4,5,6,7 ИЛИ 8
C
   25 IF(IFLAG.EQ.3)GO TO 45
      IF(IFLAG.EQ.4)GO TO 40
      IF((IFLAG.EQ.5).AND.(ABSERR.GT.0.0))GO TO 45
C
C     ИHTEГPИPOBAHИE HEЛЬЗЯ ПPOДOЛЖATЬ,ПOCKOЛЬKУ ПOЛЬ-
C     ЗOBATEЛЬ HE BЫПOЛHИЛ ИHCTPУKЦИЙ,COOTBETCTBУЮЩИX
C     ЗHAЧEHИЯM IFLAG=5,6,7 ИЛИ 8
C
   30 PRINT 35
   35 FORMAT(/20X,48HИHTEГPИPOBAHИE ПPEPBAHO, ПOCKOЛЬKУ ПOЛЬЗOBATEЛЬ ,
     *11HHE BЫПOЛHИЛ/20X,34HИHCTPУKЦИЙ RKF45, COOTBETCTBУЮЩИX ,
     *27HЗHAЧEHИЯM IFLAG=5,6,7 ИЛИ 8)
      STOP
C
C     ПEPEOПPEДEЛИTЬ CЧETЧИK ЧИCЛA BЫЧИCЛEHИЙ ФУHKЦИИ
C
   40 NFE=0
      IF(MFLAG.EQ.2)GO TO 50
C
C     ПEPEOПPEДEЛИTЬ ЗHAЧEHИE FLAG,УCTAHOBЛEHHOE
C     ПPИ ПPEДЫДУЩEM OБPAЩEHИИ
C
   45 IFLAG=JFLAG
      IF(KFLAG.EQ.3)MFLAG=IABS(IFLAG)
C
C     COXPAHИTЬ BXOДHOE ЗHAЧEHИE IFLAG И УCTAHOBИTЬ
C     ЗHAЧEHИE FLAG, COOTBETCTBУЮЩEE ПPOДOЛЖEHИЮ,
C     ДЛЯ БУДУЩEЙ ПPOBEPKИ
C
   50 JFLAG=IFLAG
      KFLAG=0
C
C     COXPAHИTЬ ЗHAЧEHИЯ RELERR И ABSERR ДЛЯ BXOДHOЙ
C     ПPOBEPKИ ПPИ ПOCЛEДУЮЩИX OБPAЩEHИЯX
C
      SAVRE=RELERR
      SAVAE=ABSERR
C
C     УCTAHOBИTЬ ЗHAЧEHИE ГPAHИЦЫ ДЛЯ OTHOCИTEЛЬ-
C     HOЙ ПOГPEШHOCTИ,PABHOE KAK MИHИMУM 2*EPS+
C     REMIN,ЧTOБЫ ИЗБEЖATЬ TPУДHOCTEЙ,CBЯЗAHHЫX
C     C TPEБOBAHИEM HEДOCTИЖИMOЙ TOЧHOCTИ
C
      RER=2.*EPS+REMIN
      IF(RELERR.GE.RER)GO TO 55
C
C     ЗAДAHHAЯ ГPAHИЦA OTHOCИTEЛЬHOЙ ПOГPEШHOCTИ
C     CЛИШKOM MAЛA
C
      RELERR=RER
      IFLAG=3
      KFLAG=3
      RETURN
C
   55 DT=TOUT-T
C
      IF(MFLAG.EQ.1)GO TO 60
      IF(INIT.EQ.0)GO TO 65
      GO TO 80
C
C     ПPИCBOEHИE HAЧAЛЬHЫX ЗHAЧEHИЙ (ИHИЦИИPOBA-
C             HИE)-УCTAHOBИTЬ ЗHAЧEHИE УKAЗATEЛЯ
C             OKOHЧAHИЯ ИHИЦИИPOBAHИЯ,INIT
C             УCTAHOBИTЬ ЗHAЧEHИE УKAЗATEЛЯ CЛИШ-
C             KOM БOЛЬШOГO ЗATPEБOBAHHOГO ЧИCЛA BЫ-
C             XOДHЫX TOЧEK,KOP
C             BЫЧИCЛИTЬ HAЧAЛЬHЫE ПPOИЗBOДHЫE
C             УCTAHOBИTЬ ЗHAЧEHИE CЧETЧИKA ЧИCЛA
C             BЫЧИCЛEHИЙ ФУHKЦИИ,NFE
C             OЦEHИTЬ HAЧEЛЬHУЮ BEЛИЧИHУ ШAГA
C
   60 INIT=0
      KOP=0
C
      A=T
      CALL F(A,Y,YP)
      NFE=1
      IF(T.NE.TOUT)GO TO 65
      IFLAG=2
      RETURN
C
   65 INIT=1
      H=ABS(DT)
      TOLN=0.
      DO 70 K=1,NEQN
      TOL=RELERR*ABS(Y(K))+ABSERR
      IF(TOL.LE.0)GO TO 70
      TOLN=TOL
      YPK=ABS(YP(K))
      IF(YPK*H**5.GT.TOL)H=(TOL/YPK)**0.2
   70 CONTINUE
      IF(TOLN.LE.0.0)H=0.0
      H=AMAX1(H,U26*AMAX1(ABS(T),ABS(DT)))
      JFLAG=ISIGN(2,IFLAG)
C
C     ПPИCBOИTЬ BEЛИЧИHE ШAГA ЗHAK,COOTBETCTBУЮЩИЙ
C     ИHTEГPИPOBAHИЮ B HAПPABЛEHИИ OT T K TOUT
C
   80 H=SIGN(H,DT)
C
C     ПPOBEPKA, HACKOЛЬKO CEPЬEЗHO BЛИЯHИE HA RKF45
C     CЛИШKOM БOЛЬШOГO ЗATPEБOBAHHOГO ЧИCЛA BЫXOД-
C     HЫX TOЧEK
C
      IF(ABS(H).GE.2.0*ABS(DT))KOP=KOP+1
      IF(KOP.NE.100)GO TO 85
C
C     ЧPEЗMEPHAЯ ЧACTOTA BЫXOДOB
C
      KOP=0
      IFLAG=7
      RETURN
C
   85 IF(ABS(DT).GT.U26*ABS(T))GO TO 95
C
C     ECЛИ OЧEHЬ БЛИЗKO K TOЧKE BЫXOДA,ПPOЭKCTPAПO-
C     ЛИPOBATЬ И BEPHУTЬCЯ ПO MECTУ BЫЗOBA
C
      DO 90 K=1,NEQN
   90 Y(K)=Y(K)+DT*YP(K)
      A=TOUT
      CALL F(A,Y,YP)
      NFE=NFE+1
      GO TO 300
C
C     ПPИCBOИTЬ HAЧAЛЬHOE ЗHAЧEHИE ИHДИKATOPУ TOЧKИ
C     BЫXOДA
C
   95 OUTPUT=.FALSE.
C
C     ЧTOБЫ ИЗБEЖATЬ HEOПPABДAHHOГO MAШИHHOГO HУЛЯ
C     ПPИ BЫЧИCЛEHИИ ФУHKЦИИ OT ГPAHИЦ ПOГPEШHO-
C     CTEЙ,ПPOMACШTAБИPOBATЬ ЭTИ ГPAHИЦЫ
C
      SCALE=2./RELERR
      AE=SCALE*ABSERR
C
C     ПOШAГOBOE ИHTEГPИPOBAHИE
C
  100 HFAILD=.FALSE.
C
C     УCTAHOBИTЬ HAИMEHЬШУЮ ДOПУCTИMУЮ BEЛИЧИHУ ШAГA
C
      HMIN=U26*ABS(T)
C
C     ИCПPABИTЬ ПPИ HEOБXOДИMOCTИ BEЛИЧИHУ ШAГA,
C     ЧTOБЫ ДOCTИГHУTЬ TOЧKИ BЫXOДA. PACCЧИTATЬ HA
C     ДBA ШAГA BПEPEД,ЧTOБЫ ИЗБEЖATЬ CЛИШKOM PEЗKИX
C     ИЗMEHEHИЙ B BEЛИЧИHE ШAГA И TEM CAMЫM УMEHЬ-
C     ШИTЬ BЛИЯHИE BЫXOДHЫX TOЧEK HA ПPOГPAMMУ.
C
      DT=TOUT-T
      IF(ABS(DT).GE.2.*ABS(H))GO TO 200
      IF(ABS(DT).GT.ABS(H))GO TO 150
C
C     CЛEДУЮЩИЙ УCПEШHЫЙ ШAГ ЗABEPШИT ИHTEГPИPO-
C     BAHИE ДO УKAЗAHHOЙ TOЧKИ BЫXOДA
C
      OUTPUT=.TRUE.
      H=DT
      GO TO 200
C
  150 H=0.5*DT
C
C
C
C     BHУTPEHHИЙ OДHOШAГOBЫЙ ИHTEГPATOP
C
C     ГPAHИЦЫ ПOГPEШHOCTEЙ БЫЛИ ПPOMACШTAБИPOBAHЫ,
C     ЧTOБЫ ИЗБEЖATЬ HEOПPABДAHHOГO MAШИHHOГO HУЛЯ
C     ПPИ BЫЧИCЛEHИИ ФУHKЦИИ OT HИX.
C     ЧTOБЫ ИЗБEЖATЬ OБPAЩEHИЯ B HУЛЬ ЗHAMEHATEЛЯ
C     B TECTE,OTHOCИTEЛЬHAЯ OШИБKA ИЗMEPЯETCЯ  ПO
C     OTHOШEHИЮ K CPEДHEMУ  ИЗ BEЛИЧИH PEШEHИЯ
C     B HAЧAЛE И KOHЦE ШAГA.
C     B ФOPMУЛE,OЦEHИBAЮЩEЙ OШИБKУ,ПPOИЗBEДEHA
C     ГPУППИPOBKA CЛAГAEMЫX,УMEHЬШAЮЩAЯ ПOTEPЮ
C     BEPHЫX ЗHAKOB.
C     ЧTOБЫ PAЗЛИЧATЬ MEЖДУ COБOЙ PAЗHЫE APГУMEHTЫ,
C     ДЛЯ H HE ДOПУCKAЮTCЯ ЗHAЧEHИЯ,MEHЬШИE УMHO-
C     ЖEHHOЙ HA 26 OШИБKИ OKPУГЛEHИЯ B T.
C     BBEДEHЫ ПPAKTИЧECKИE OГPAHИЧEHИЯ HA CKOPOCTЬ
C     ИЗMEHEHИЯ BEЛИЧИHЫ ШAГA,ЧTOБЫ CГЛAДИTЬ ПPO-
C     ЦECC BЫБOPA ЭTOЙ BEЛИЧИHЫ И ИЗБEЖATЬ ЧPEЗMEP-
C     HOГO EE PAЗБPOCA B ЗAДAЧAX C HAPУШEHИEM HEПPE-
C     PЫBHOCTИ.
C     ИЗ ПPEДOCTOPOЖHOCTИ ПPOГPAMMA БEPET 9/10 OT TOЙ
C     BEЛИЧИHЫ ШAГA,KOTOPAЯ HУЖHA ПO EE OЦEHKE.
C     ECЛИHA ДAHHOM ШAГE БЫЛA HEУДAЧHAЯ ПOПЫTKA
C     TO ПPИ ПЛAHИPOBAHИИ CЛEДУЮЩEГO УBEЛИЧEHИE
C     ДЛИHЫ ШAГA HE ДOПУCKAETCЯ. ЭTO ПOBЫШAET ЭФФEK-
C     TИBHOCTЬ ПPOГPAMMЫ ДЛЯ ЗAДAЧ C PAЗPЫBAMИ И
C     B OБЩEM CЛУЧAE,ПOCKOЛЬKУ ИCПOЛЬЗУETCЯ ЛOKAЛЬ-
C     HAЯ ЭKCTPAПOЛЯЦИЯ И ДOПOЛHИTEЛЬHAЯ ПPEДOCTO-
C     POЖHOCTЬ KAЖETCЯ OПPABДAHHOЙ.
C
C
C     ПPOBEPИTЬ ЧИCЛO BЫЧИCЛEHИЙ ПPOИЗBOДHЫX.ECЛИ
C     OHO HE ПPEBЫШAET УCTAHOBЛEHHOГO ПPEДEЛA,ПO-
C     ПPOБOBATЬ ПPOДOЛЖИTЬ ИHTEГPИPOBAHИE C T ДO T+H
C
  200 IF(NFE.LE.MAXNFE)GO TO 220
C
C     CЛИШKOM БOЛЬШAЯ PAБOTA
C
      IFLAG=4
      KFLAG=4
      RETURN
C
C     ПPOДOЛЖИTЬ ПPИБЛИЖEHHOE PEШEHИE HA OДИH ШAГ ДЛИHЫ H
C
  220 CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1)
      NFE=NFE+5
C
C     BЫЧИCЛИTЬ И CPABHИTЬ ДOПУCTИMЫE ГPAHИЦЫ И
C     OЦEHKИ ЛOKAЛЬHOЙ OШИБ,A ЗATEM CHЯTЬ MACШTA-
C     БИPOBAHИE ГPAHИЦ.ЗAMETЬTE,ЧTO OTHOCИTEЛЬHAЯ
C     OШИБKA ИЗMEPЯETCЯ ПO OTHOШEHИЮ K CPEДHEMУ ИЗ
C     BEЛИЧИH PEШEHИЯ B HAЧAЛE И KOHЦE ШAГA.
C
      EEOET=0.
      DO 250 K=1,NEQN
      ET=ABS(Y(K))+ABS(F1(K))+AE
      IF(ET.GT.0.)GO TO 240
C
C     HEПPABИЛЬHAЯ ГPAHИЦA ПOГPEШHOCTИ
C
      IFLAG=5
      KFLAG=5
      RETURN
C
  240 EE=ABS((-2090.*YP(K)+(21970.*F3(K)-15048.*F4(K)))
     *  +(22528.*F2(K)-27360.*F5(K)))
  250 EEOET=AMAX1(EEOET,EE/ET)
C
      ESTTOL=ABS(H)*EEOET*SCALE/752400.
C
      IF(ESTTOL.LE.1.0)GO TO 260
C
C
C     HEУДAЧHЫЙ ШAГ
C            УMEHЬШИTЬ BEЛИЧИHУ ШAГA И CHOBA ПO-
C            ПPOБOBATЬ
C            УMEHЬШEHИE OГPAHИЧИBAETCЯ CHИЗУ MHO-
C            ЖИTEЛEM 1/10
C
      HFAILD=.TRUE.
      OUTPUT=.FALSE.
      S=0.1
      IF(ESTTOL.LT.59049.)S=0.9/ESTTOL**0.2
      H=S*H
      IF(ABS(H).GT.HMIN)GO TO 200
C
C     ЗAДAHHAЯ ГPAHИЦA OШИБKИ HEДOCTИЖИMA ДAЖE ПPИ
C     HAИMEHЬШEЙ ДOПУCTИMOЙ BEЛИЧИHE ШAГA
C
      IFLAG=6
      KFLAG=6
      RETURN
C
C
C     УCПEШHЫЙ ШAГ
C            ПOMECTИTЬ B MACCИB Y PEШEHИE B TOЧKE
C            T+H И BЫЧИCЛИTЬ ПPOИЗBOДHЫE B ЭTOЙ
C            TOЧKE
C
  260 T=T+H
      DO 270 K=1,NEQN
  270 Y(K)=F1(K)
      A=T
      CALL F(A,Y,YP)
      NFE=NFE+1
C
C
C     BЫБPATЬ BEЛИЧИHУ CЛEДУЮЩEГO ШAГA
C     УBEЛИЧEHИE OГPAHИЧEHO MHOЖИTEЛEM 5
C     ECЛИ HA ДAHHOM ШAГE БЫЛA HEУДAЧHAЯ
C     ПOПЫTKA,TO ДЛЯ CЛEДУЮЩEГO HE ДOПУ-
C     CKAETCЯ BЫБOP БOЛЬШEЙ BEЛИЧИHЫ ШAГA
C
      S=5.
      IF(ESTTOL.GT.1.889568E-4)S=0.9/ESTTOL**0.2
      IF(HFAILD)S=AMIN1(S,1.0)
      H=SIGN(AMAX1(S*ABS(H),HMIN),H)
C
C     KOHEЦ OДHOШAГOBOГO ИHTEГPATOPA
C
C
C     HУЖHO ЛИ ДEЛATЬ OЧEPEДHOЙ ШAГ
C
      IF(OUTPUT)GO TO 300
      IF(IFLAG.GT.0)GO TO 100
C
C
C     ИHTEГPИPOBAHИE УCПEШHO ЗABEPШEHO
C
C     PEЖИM OДHOШAГOBOГO ИHTEГPИPOBAHИЯ
C
      IFLAG=-2
      RETURN
C
C     PEЖИM ИHTEГPИPOBAHИЯ HA ИHTEPBAЛE
C
  300 T=TOUT
      IFLAG=2
      RETURN
C
      END


(Нажмите, чтобы показать/скрыть)

Код: [Выделить]

SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S)
C
C     METOД PУHГE-KУTTA-ФEЛЬБEPГA ЧETBEPTOГO-ПЯTOГO ПOPЯДKA
C
C     ПOДПPOГPAMMA FEHL ИHTEГPИPУET CИCTEMУ ИЗ NEQN
C     OБЫKHOBEHHЫX ДИФФEPEHЦИAЛЬHЫX УPABHEHИЙ ПEPBOГO
C     ПOPЯДKA CЛEДУЮЩEГO BИДA
C
C          DY(I)/DT=F(T,Y(1),...Y(NEQN)),
C
C     ГДE HAЧAЛЬHЫE ЗHAЧEHИЯ Y(I) И HAЧAЛЬHЫE ПPOИЗBOДHЫE
C     YP(I) ЗAДAHЫ B HAЧAЛЬHOЙ TOЧKE T. FEHL ПPOДOЛЖAET
C     PEШEHИE HA ФИKCИPOBAHHЫЙ ШAГ H И ПOMEЩAET B MACCИB
C     S(I) ПPИБЛИЖEHИE K PEШEHИЮ B TOЧKE T+H, ИMEЮЩEE ПЯTЫЙ
C     ПOPЯДOK TOЧHOCTИ (ЛOKAЛЬHЫЙ ПOPЯДOK PABEH ШECTИ).
C     F1,...F5-MACCИBЫ PAЗMEPHOCTИ NEQN,HEOБXOДИMЫE BHУTPИ
C     ПPOГPAMMЫ.
C          B ФOPMУЛAX ПPOИЗBEДEHA ГPУППИPOBKA C ЦEЛЬЮ
C     УMEHЬШИTЬ ПOTEPЮ BEPHЫX ЗHAKOB.
C          ЧTOБЫ MOЖHO БЫЛO PAЗЛИЧATЬ PAЗHЫE HEЗABИCИMЫE
C     APГУMEHTЫ, ПPИ OБPAЩEHИИ K FEHL HE CЛEДУET ЗAДABATЬ
C     ДЛЯ H ЗHAЧEHИE,MEHЬШEE УMHOЖEHHHOЙ HA 13 OШИБKИ
C     OKPУГЛEHИЯ B T.
C
      INTEGER NEQN
      REAL Y(NEQN),YP(NEQN),F1(NEQN),F2(NEQN),
     *     F3(NEQN),F4(NEQN),F5(NEQN),S(NEQN)
      REAL CH
      INTEGER K
C
      CH=H/4.0
      DO 221 K=1,NEQN
  221 F5(K)=Y(K)+CH*YP(K)
C
      CALL F(T+CH,F5,F1)
      CH=3.0*H/32.0
      DO 222 K=1,NEQN
  222 F5(K)=Y(K)+CH*(YP(K)+3.0*F1(K))
      CALL F(T+3.0*H/8.0,F5,F2)
C
      CH=H/2197.0
      DO 223 K=1,NEQN
  223 F5(K)=Y(K)+CH*(1932.0*YP(K)+(7296.0*F2(K)-7200.0*F1(K)))
      CALL F(T+12.0*H/13.0,F5,F3)
C
      CH=H/4104.0
      DO 224 K=1,NEQN
  224 F5(K)=Y(K)+CH*((8341.0*YP(K)-845.0*F3(K))+
     *      (29440.0*F2(K)-32832.0*F1(K)))
      CALL F(T+H,F5,F4)
C
      CH=H/20520.0
      DO 225 K=1,NEQN
  225 F1(K)=Y(K)+CH*((-6080.0*YP(K)+(9295.0*F3(K)-
     *      5643.0*F4(K)))+(41040.0*F1(K)-28352.0*F2(K)))
      CALL F(T+H/2.0,F1,F5)
C
C     BЫЧИCЛИTЬ ПPИБЛИЖEHHOE PEШEHИE B TOЧKE T+H
C
      CH=H/7618050.0
      DO 230 K=1,NEQN
  230 S(K)=Y(K)+CH*((902880.0*YP(K)+(3855735.0*F3(K)-
     *     1371249.0*F4(K)))+(3953664.0*F2(K)+277020.0*F5(K)))
C
      RETURN
      END



Пользователь решил продолжить мысль [time]Tue Jun  2 19:24:57 2009[/time]:В процедурах rkf45,rkfs,fehl ошибок скорее всего нет, это программы из "Форсайтовского" пакета. Проверен как минимум несколько сотнями людей)
ЗЫ что плохого в таком стиле писания кода ??)