Методы решения обычных дифференциальных уравнений (Лабораторная работа № 7), страница 6

C        Запись информацыи для выбраннога узла в выходной файл

С       AAAA.REZ

WRITE(3,103)X,(Y(I),I=1,N)

C        Запись информации для избранного узла в выходной файл

С       AAAA.GRF

WRITE(4,103)X,(Y(I),I=1,N)

103 FORMAT(10F8.4)

C        Вывод информации  на экран дисплея для текущего контроля

WRITE(5,106)X,NREZ

106 FORMAT(1X,'X=',F8.4,'    NREZ=',I4)

C        Обращение к подпрограмме решения дифференциальных уравнений

6 CALL EILER(N,X,H,Y,PRAV)

C        Надо ли  продолжать расчет?

9 IF(NREZ.LT.501)GO TO 10

WRITE(5,*)'DIFEIL: РАСЧЕТ ПРИОСТАНОВЛЕН ПО ПРИЧИНЕ',

*' ЗАПОЛНЕНИЯ ВЫХАДНЫХ ФАЙЛАУ (NREZ>501)'

GO TO 11

10 IF ((X-XKAN).GT.0.)GO TO 11

C        Увеличение независимой переменной на величину шага H

X=X+H

GO TO 4

C        Закрытие выходных файлов

11 CLOSE(UNIT=3)

CLOSE(UNIT=4)

WRITE(5,101) Z2,Z3               

101 FORMAT(1X,'DIFEIL: РАБОТА ЗАКОНЧЕНА.'

*' Р-ТЫ ЗАПИСАНЫ В ФАЙЛЫ ',A8,' I ',A8)

STOP

END

C

SUBROUTINE EILER(N,X,H,Y,FNC)

C ПОДПРОГРАММА РЕШЕНИЯ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ

С УРАВНЕНИЙ dY/dx=F(x,Y)

C С НАЧАЛЬНЫМИ УСЛОВИЯМИ (ЗАДАЧА КАШИ)

C ЯВНЫМ МЕТОДОМ ЭЙЛЕРА НА ОДНОМ ВРЕМЕННОМ ШАГЕ

C ОБОЗНАЧЕНИЯ ПЕРЕМЕННЫХ:

C    N - КОЛИЧЕСТВО УРАВНЕНИЙ В СИСТЕМЕ (ВХОДНАЯ

С  ВЕЛИЧИНА)

C    X - ТЕКУЩЕЕ ЗНАЧЕНИЕ НЕЗАВИСИМОЙ ПЕРЕМЕННОЙ x

C        (РАБОЧАЯ И ВЫХАДНАЯ ВЕЛИЧИНА)

C    H - ШАГ ИЗМЕНЕНИЯ НЕЗАВИСИМОЙ ПЕРЕМЕННОЙ x

C        (ВХОДНАЯ ВЕЛИЧИНА)

C    Y(N) - МАССИВ ТЕКУЩИХ ЗНАЧЕНИЙ  ИНТЕГРИРУЕММЫХ

С   ПЕРЕМЕННЫХ Y

C        (РАБОЧИЙ МАССИВ)

C    F(N) - МАССИВ ЗНАЧЕНИЙ ФУНКЦИЙ F(x,Y) В ПРАВЫХ ЧАСТЯХ

C УРАВНЕНИЙ (РАБОЧИЙ МАССИВ).ЗНАЧЕНИЯ ФУНКЦИЙ

С ВЫССЧИТЫВАЮТСЯ  ВО ВНЕШНЕЙ ПОДПРОГРАММЕ FNC,

C КОТОРАЯ ЦЕЛИКОМ ОПРЕДЕЛЯЕТСЯ ВИДОМ

C ПРАВОЙ ЧАСТИ СИСТЕМЫ УРАВНЕНИЙI ДЛЯ КАЖДОЙ НОВОЙ

C СИСТЕМЫ  СКЛАДЫВАЕТСЯ ПОЛЬЗОВАТЕЛЕМ ЗАНОВО 

C     FNC - ИМЯ ВНЕШНЕЙ ПОДПРОГРАММЫ, В КОТОРОЙ

С ВЫСЧИТЫВАЮТСЯ ЗНАЧЕНИЯ ФУНКЦИЙ F(x,Y) В ПРАВОЙ

C ЧАСТИ СIСТЭМЫ.

C ИМЯ ЭТОЙ ПОДПРОГРАММЫ В СПИСКЕ ФОРМАЛЬНЫХ

C ПАРАМЕТРОВ ОПЕРАТОРА SUBROUTINE ЭТОЙ ПОДПРОГРАММЫ

C И В СПИСКЕ ФАКТИЧЕСКИХ  ПАРАМЕТРОВ ОПЕРАТОРА CALL

С ПРОГРАММА, КОТОРАЯ ЕЕ ВЫЗЫВАЕТ,

C МОЖЕТ БЫТЬ РАЗНОЙ, В ТАКОМ СЛУЧАЕ В ПРОГРАММЕ,

С КОТОРАЯ ЕЕ ВЫЗЫВАЕТ,

C ЭТО ИМЯ ДОЛЖНО БЫТЬ ОПИСАНО В ОПЕРАТОРЕ EXTERNAL

DIMENSION Y(*),F(10)

C ВЫЧИСЛЕНИЕ ЗНАЧЭНИЙ ФУНКЦИЙ F(x,Y) В ПРАВЫХ ЧАСТЯХ

С УРАВНЕНИЙ

DO 5 I=1,N

5 CALL FNC(X,Y,F)

C ВЫЧЕСЛЕНИЕ ЗНАЧЕНИЙ  ИНТЕГРИРУЕМЫХ ПЕРЕМЕННЫХ Y

С НА НОВОМ ШАГЕ

DO 6 I=1,N

6 Y(I)=Y(I)+H*F(I)

11 RETURN

END

C

C ПОДПРОГРАММА ВЫЧИСЛЕНИЯ ПРАВЫХ ЧАСТЕЙ СИСТЕМЫ

С ДЛЯ СЛЕДУЮЩЕЙ СIСТЭМЫ ДИФФЕРЕНЦИАЛЬНЫХ

C УРАВНЕНИЙ 3-ГО ПОРЯДКА:

C   dy1/dx=f1(y1,y2,y3,x); где: f1(y1,y2,y3,x)=cos(x)-exp(-x)-y3

C   dy2/dx=f2(y1,y2,y3,x);      f2(y1,y2,y3,x)=y3

C   dy3/dx=f3(y1,y2,y3,x);      f3(y1,y2,y3,x)=-y2

SUBROUTINE PRAV(X,Y,F)

DIMENSION Y(*),F(*)

F(1)=COS(X)-EXP(-X)-Y(3)

F(2)=Y(3)

F(3)=-Y(2)

RETURN

END

Дополнение 2

C  ПРОГРАММА 'DIFRK2' 16.09.94

C  РЕШЕНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

С dY/dx = F(x,Y),

C   Y(0)=Y0 МЕТОДОМ РУНГЕ-КУТТА ВТОРОГО ПОРЯДКА ПО  

С   АЛГОРИТМУ:

C         y(i,j+1)=y(i,j)+(h/2)*(k1(i,j)+k2(i,j)),

C         k1(i,j)=f(i)(x(j),y(i,j)),

C         k2(i,j)=f(i)(x(j)+h,y(i,j)+h*k1(i,j)),

C         i=1,2, ... ,n (n-КОЛИЧЕСТВО УРАВНЕНИЙ В СИСТЕМЕ)

C         j=0,1,2, ... ,m-1 (m-КОЛИЧЕСТВО УЗЛОВ ПО ОСИ x)

C      РАЗРАБОТАНА НА КАФЕДРЕ ЭЛЕКТРИЧЕСКИХ СТАНЦИЙ

С    20.09.94.

C            ОБОЗНАЧЕНИЯ ПЕРЕМЕННЫХ В ПРОГРАММЕ:

C N - КОЛИЧЕСТВО ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В

С  СИСТЕМЕ (ВХОДНАЯ ВЕЛИЧИЧНА). МАКСИМАЛЬНОЕ

C   КОЛИЧЕСТВО УРАВНЕНИЙ СИСТЕМЫ ПРИНЯТО

C   В ПРОГРАММЕ РАВНОЕ 16, НО ЭТО ОГРАНИЧЕНИЕ СВЯЗАНО

С ТОЛЬКО С ПРИНЯТЫМИ В ПРОГРАММЕ

C   РАЗМЕРНОСТЯМИ МАССИВОВ Y(0), Y0(0) И

C   РАЗМЕРАМИ ТАБЛИЦ РЕЗУЛЬТАТОВ AAAA.REZ I AAAA.GRF

C INTWR - ИНТЕРВАЛ ПЕЧАТИ-КОЛИЧЕСТВО ШАГОВ ПО

С НЕЗАВИСИМОЙ ПЕРЕМЕННОЙ x, КОТОРЫЕ ПРОПУСКАЮТСЯ

C   ПРИ ВЫВОДЕ РЕЗУЛЬТАТОВ РЕШЕНИЯ НА ЭКРАН ИЛИ

C   ПРИНТЕР (ВХОДНАЯ ВЕЛИЧИНА). НАПРИМЕР, ПРИ X0=0, XKAN=0.1,

C   H=0.0001 РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ  БУДУТ

С ВЫСЧИТЫВАТЬСЯ В (XKAN- X0)/H=(0.1-0.0)/0.0001=1000