4.6 Блок - схема алгоритма метода Рунге - Кутта с автоматическим выбором шага

Рисунок 3

4.7 Подпрограмма метода Рунге - Кутта с автоматическим выбором шага

Подпрограмма ARK позволяет решать произвольную систему N-го порядка с автоматическим выбором шага интегрирования. Эта подпрограмма обращается:

· к подпрограмме одного шага- SH,

· к подпрограмме вычисления правых частей системы,

· к подпрограмме вывода.

Подпрограмма SH записана в универсальном виде и приведена выше. Головной вызывающий модуль, а также подпрограммы правых частей и вывода Пользователь должен составить самостоятельно.

В главном модуле оператором EXTERNAL должны быть объявлены имена подпрограмм правых частей и вывода. Оператором DIMENSION должны быть объявлены массивы - фактические параметры подпрограммы ARK. Эти массивы, по желанию, могут объявляться как одномерные или как двухмерные. Размеры массивов (N,4),(N,3),(N,4), где N-порядок системы. Формальные имена этих массивов в подпрограмме ARK, соответственно, X,R,F. В главном модуле первые N элементов массива, соответствующего X, заполняются начальными условиями, а следующие N элементов заполняются весовыми коэффициентами погрешности. В подпрограммах правых частей и вывода в первых N элементах массива, соответствующего X, при входе содержатся текущие значения всех N переменных системы и не должны там переопределяться. Первые N элементов массива, соответствующего F, должны заполняться в подпрограмме правых частей вычисляемыми там значениями правых частей системы.

Формальными параметрами в подпрограмме правых частей должны быть параметры (T,X,F,N), где T-независимая переменная системы.

Формальными параметрами подпрограммы вывода должны быть параметры (T,X,F,N,IER), где IER- код ошибки, определяемый в подпрограмме ARK:

· IER=0,- ошибки нет;

· IER=1,- знак заданного начального шага не соответствует движению от начала интервала интегрирования к его концу;

· IER=2,- начальный шаг или/и длина интервала интегрирования ошибочно заданы равными нулю;

· IER=3,- шаг в процессе счёта стал более чем в 1000 раз меньше начального.

Массивы X и F в подпрограммах правых частей и вывода можно объявлять как одномерные, с регулируемым размером X(N),F(N).

В главном модуле для подпрограммы ARK должны задаваться максимальный (он же и начальный) шаг интегрирования HM, начало TN и конец TK интервала интегрирования, а также значение требуемой абсолютной погрешности решения E.

Подпрограмма ARK вычисляет решение системы и в каждой точке, удовлетворяющей условиям точности, обращается к подпрограмме вывода, передавая ей значения параметров T,X,F,IER. Пользователь может запрограммировать здесь печать необходимых переменных или накопление их в дополнительных массивах для последующей обработки. (В последнем случае дополнительные массивы следует передавать в главный модуль через общую область памяти с помощью оператора COMMON). После возврата из подпрограммы вывода, ARK продолжает вычисление следующей точки решения.

SUBROUTINE ARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)

C Подпрограмма автоматического выбора шага.

C HM -Задаваемый максимальный шаг.

C TN,TK -Начало и конец отрезка интегрирования.

C N -Порядок системы.

C E -Задаваемое значение абсолютной погрешности.

EXTERNAL PRAV,OUT

C PRAV и OUT имена составляемых Пользователем подпрограмм правых частей и вывода.

C IER -Выходной код ошибки.

DIMENSION X(N,4),R(N,3),F(N,4)

C Первый столбец массива X должен при входе содержать начальные условия,

С на выходе в нем содержится решение.

C Второй столбец массива X должен при входе содержать весовые коэффициенты погрешности.

C Первый столбец массива F должен заполняться вычисляемыми

C в подпрограмме PRAV значениями правых частей системы уравнений.

C Остальные элементы массивов X,R,F -рабочие.

DO 3 K=1,N

3 F(K,4)=X(K,2)

T=TN

HB=2*HM

IER=0

KP=0

KLP=1

CALL PRAV(T,X,F,N)

C Вызов составленной Пользователем подпрограммы правых частей системы уравнений.

C T -Независимая переменная системы.

IF((TK-TN)*HM)4,5,60

4 IER=1

GO TO 60

5 IER=2

60 CALL OUT(T,X,F,N,IER)

C Вызов составленной Пользователем подпрограммы вывода результатов шага

IF(IER.NE.0)RETURN

6 H=HB/2

CALL SH(T,H,X,X(1,2),F(1,2),PRAV,N,R)

8 T1=T+H

CALL SH(T1,H,X(1,2),X(1,3),F(1,3),PRAV,N,R)

T2=T+HB

CALL SH(T,HB,X,X(1,4),F,PRAV,N,R)

EH=0

DO 14 K=1,N

14 EH=EH+ABS((X(K,3)-X(K,4))/15*F(K,4))

IF(EH-E)21,21,16

16 HB=HB/2

IF(HB.LT.HM/512)IER=3

IF(IER.EQ.3)RETURN

KP=0

GO TO 6

21 T=T1

IF(KLP.EQ.1)CALL OUT(T1,X(1,2),F(1,2),N,IER)

KLP=0

CALL OUT(T2,X(1,3),F(1,3),N,IER)

IF(KP.EQ.1)RETURN

IF(ABS(TK-T2)-ABS(HB))29,29,30

29 KP=1

HB=TK-T2

GO TO 33

30 IF(EH-E/50)31,37,37

31 IF(2*H.LE.-HM.AND.ABS(TK-T2).GE.ABS(2*HB))GO TO 32

37 DO 38 K=1,N

X(K,1)=X(K,2)

38 X(K,2)=X(K,3)

GO TO 8

32 HB=2*HB

33 T=T2

DO 35 K=1,N

35 X(K,1)=X(K,3)

KLP=1

GO TO 6

END

4.8 Тестовая задача

Решим дифференциальное уравнение

с начальными условиями . Легко видеть, что решением этой задачи является функция . Вычислим решение на интервале , что составит почти периода этой функции. Если при таком длительном интегрировании амплитуда косинусоиды существенно не изменится, то алгоритм численно устойчив. Можно также сравнить решение в конечной точке

Подпрограмма правых частей для этого уравнения будет такой.

SUBROUTINE PRAV(T,X,F,N)

DIMENSION X(N),F(N)

F(1)=X(2)

F(2)= -X(1)

RETURN

END

В подпрограмме вывода предусмотрим заполнение результатами массива D для построения графиков на интервале в пять периодов, а также заполнение массива С положительными максимумами вычисляемой функции на всем интервале интегрирования. Эти массивы передадим в главный модуль через общую область.

SUBROUTINE OUT(T,X,F,N,IER)

DIMENSION X(N),F(N),D(3,1000),C(300)

COMMON K,L,KP,D,C

IF(T.LT.31.4)THEN

K=K+1