Теоретическая часть
Программа предназначена для численного решения системы
обыкновенных дифференциальных уравнений вида:
Y'=F(X,Y), с начальными условиями Y(X0)=Yo на отрезке
[X,X] методом Хемминга с постоянным шагом интегрирования. В
каждой i+1 точке находим начальное приближение Р к решению Y
по предсказывающей формуле:
Pi+1=Yi-3+4*h*(2*Y'i-Y'i-1+2*Y'i-2)/3, где
Yi-3 решение в i-3 точке,
Y'i,Y'i-1,Y'i-2 - значения производных в точках i, i-1,
i-2 соответственно.
Для улучшения решения используется корректирующая формула
Ci+1=[9*Yi-Yi-2+3*h*(M'i+1+2*Y'i-Y'i-1)]/8, где
Mi+1=Pi+1-112*(P-Ci)/121; M'i+1=F(Xi+1,Mi+1).
Решение системы в i+1 точке находится по формуле
G=Wj*|Pi+1,j-Ci+1|, где
Wj=1
j- номер компоненты вектора.
На участке "разгона" значения Yi-k и Y'i-k (k=0, 1, 2)
вычисляются методом Рунге-Кутта по формуле
Yi=Ui(2)-(Ui(i)-Ui(2))/15, где i- номер точки, в которой
ищется решение, Ui- решение системы в i-ой точке, полученное с
шагом h/l;
U(l)i-m+1/l=A(l)i-m+(K(l)1+2*K(l)2+2*K(l)3+K(l)4)i--m+1/l/6,
где
j=1, 2, ..., n,
K(l)1=h*F(Xi-m,A(l)i-m)/l;
K(l)2=h*F(Xi-m+h/2*l,A(l)i-m+K(l)1/2)/l;
K(l)3=h*F(Xi-m+h/2*l,A(l)i-m+K(l)2/2)/l;
K(l)4=h*F(Xi-m+h/l,A(l)i-m+K(l)3)/l.
A, U ,K - векторы n-го порядка; l=1, 2; m=1 при l=1; m=1,
1/2 при l=2;
A(l)i-1=Y(l)i-1; A(2)i-1/2=U(2)i-1/2.
Характеристика программы.
Программа состоит из стандартной информативы, реализующей
описанный метод, рабочей информативы, задающей правые части
уравнений системы и директивы.
Длина стандартной информативы 1600 символов. Объем исход-
ных данных : 7 чисел, 2 массива, n функций. В результате рабо-
ты программы на печать выводится на участке "разгона" X, зна-
чения функций и производных, далее X, G и Y[n] на всем отрезке
интегрирования через Ю шагов и в конце отрезка.
Программа рекомендуется для решения систем обыкновенных
дифференциальных уравнений на больших отрезках, так как счита-
ет быстрее одноточечных методов. Для контроля постоянно выво-
дится погрешность вычислений G, которая позволяет следить за
точностью решения.
"Разгон" (нахождение значений функций и производных в
точках X0, X0+Q, X0+2*Q , X0+3*Q, где Q - шаг интегрирования )
осуществляется методом Рунге-Кутта с увеличенной разрядностью.
В программе предусмотрена возможность при получении боль-
шой погрешности вычисления в точка "разгона" уменьшить шаг ин-
тегрирования в этих точках (см. способ задания J), а при быст-
ром возрастании погрешности вычислений G уменьшить шаг интег-
рирования методом Хемминга или увеличить разрядность вычисле-
ний.
Программа позволяет производить интегрирование как с по-
ложительным, так и с отрицательным шагом (соответственно меня-
ются X0, Xк и Q).
Порядок решения задачи.
Для решения задачи вводятся стандартная и рабочая инфор-
мативы и директива.
В рабочей информативе после метки Ц программа вычисления
правых частей системы. Здесь Z[1]=...; Z[2]=...; ...;
Z[n]=...; - правые части исходной системы обыкновенных диффе-
ренциальных уравнений как функции от X1 и Y[1], Y[2], ...,
Y[n], X1 - соответствует аргументу, Y[I] - соответствует функ-
циям. I=1, 2, ..., N. Операторная часть рабочей информативы
заканчивается оператором перехода "НА" Ф.
В описательной части рабочей информативы задаются X0, XK
- соответственно начало и конец отрезка интегрирования, Q -
шаг интегрирования методом Хемминга, J - число, определяющее,
во сколько раз следует уменьшить шаг интегрирования методом
Рунге-Кутта на участке "разгона" для получения решения того же
порядка точности, что и в методе Хемминга,
N=n - порядок системы;
Y[n] - вектор начальных условий,
W[n] - вектор коэффициентов для вычисления невязки
W[I]=1, и описаны
A[n], B[n], C[n] - массивы значений функций в точках i-3,
i-2, i-1 соответственно,
Я[n], Б[n], Г[n], D[n] - массивы значений производных в
точках i-3, i-2, i-1, i соответственно,
Z[n] - массив правых частей,
П[n], P[n] - рабочие массивы.
В директиве задаются : R - разрядность вычислений по ме-
тоду Хемминга ("разгон" происходит с увеличенной разряд-
ностью), Ю - число, определяющее период печати (количество ша-
гов). Директива должна оканчиваться оператором "НА" HMG.
Описание работы программы
Данная расчетно-графическая работа (далее РГР) составлена
на языке PC MathLab ( PC-MATLAB (c) Copyright The MathWorks,
Inc. 1984-1989 Version 3.5f 17-July-89 Serial Number 22961) и
выполнена в виде двух модулей (третий - контрольный пример),
распечатка которых приведена в приложении.
1. Hemming.m
"Стандартный" головной модуль.
Входные данные: отсутствуют.
Выходные данные: отсутствуют.
Язык реализации: PC MathLab.
Операционная система: MS-DOS 3.30 or higher
Пояснения к тексту модуля:
Структура данного модуля элементарна. Вначале очищается
экран, задаются исходные данные для второго модуля, как X0,XK
- начальное и конечное значение, Q - шаг, J - число, определя-
ющее во сколько раз нужно уменьшать шаг интегрирования методом
Рунге-Кутта (далее Р-К) на участке "разгона" для получения то-
го же порядка точности, что и в методе Хемминга, N - порядок
системы, Y - вектор начальных значений, W - вектор коэффициен-
тов для вычисления невязки и т.д. Затем вызывается модуль ре-
шения системы в формате:
[x,y,dg]=hem('ours',x0,xk,q,j,n,y,w,ur), где
x,y - точки решения
dg - ошибка остальные параметры описаны выше.
Необходимо отметить, что несмотря на отсутствие входных и
выходных данных, внутри данного модуля задаются начальные зна-
чения и выводятся результаты вычислений в числовом виде и гра-
фиков, а также оценка по быстродействию (TIME) и количеству
выполненных операций (FLOPS), однако эти данные нельзя охарак-
теризовать как входные и выходные.
2. Hem.m
Модуль, которые непосредственно и решает систему ОДУ ме-
тодом Хемминга.
Входные данные:
FunFcn - имя функции, вычисляющей левые части
X0,XK - начальное и конечное значение для счета
Q - шаг интегрирования
J - число, определяющее во сколько раз нужно уменьшать
шаг интегрирования методом Рунге-Кутта (далее Р-К) на участке
"разгона" для получения того же порядка точности, что и в ме-
тоде Хемминга
N - порядок системы
Y - вектор начальных значений
W - вектор коэффициентов для вычисления невязки
UR - число, определяющее период печати
Выходные данные:
x - матрица точек, для которых вычислено решение
y - матрица решений
dg - ошибка интегрирования
Язык реализации: PC MathLab.
Операционная система: MS-DOS 3.30 or higher
Пояснения к тексту модуля:
Данный модуль содержит в своем теле всего одну функцию,
входные и выходные данные которой являются входными и выходны-
ми данными текущего модуля. Они описаны выше. Мы же займемся
описанием данной функции:
После описания функции HEM устанавливается формат выход-
ных данных LONG E, а также происходит инициализация рабочих
массивов, как массивы значений функции в точках i-3, i-2, i-1;
массивы значений производных в этих же точках, массивы правых
частей и т.д. Всвязи с отсутствием в языке MathLab конструкции
безусловного перехода, используется конструкции While 1 (бес-
конечный цикл), Break (переход к началу While) и IF (Если).
Из-за таких немного "странных" конструкций вся дальнейшая
часть программы может быть весьма условно представлена такой
схемой:
While (не конец расчетов)
While 1
...
IF
...
...
END
...
...
IF
...
...
...
END
...
end
end
В целом, в данных циклах вычисляется номер шага, и если
он меньше 5, то вычисления сводятся к вычислению методом Р-К,
а если больше 5, то производятся вычисления по методу Хемминга
со всеми своими дополнительными действиями, как вычисление по
корректирующей формуле и т.д. В обоих случаях происходит об-
новление рабочих и других промежуточных массивов и вывод ин-
формации на экран. В случае решения в точках "разгона" вычис-
ляются также коэффициенты K1, K2, K3, K4, используемые в мето-
де Р-К. Также функция сама проверяет точность вычислений и в
случае необходимости корректирует шаг. Если шаг "сделан", то
программа выводит результаты на экран и заносит их в массив,
который представлен в виде нескольких столбцов. Также в необ-
ходимых для функции случаях она обращается к подпрограмме
FunFcn, которая занимается вычислением левых частей, вызов и
возврат значений которой должен быть следующим:
Z=feval(FunFcn,x,y), где
Z - вектор вычисленных левых частей,
X,Y - векторы точек, для которых