Программа предназначена для численного решения системы
обыкновенных дифференциальных уравнений вида:
Y'=F(X,Y), с начальными условиями Y(X 40 0)=Y 4o 0 на отрезке
[X,X] методом Хемминга с постоянным шагом интегрирования. В
каждой i+1 точке находим начальное приближение Р к решению Y
по предсказывающей формуле:
P 4i+1 0=Y 4i-3 0+4*h*(2*Y' 4i 0-Y' 4i-1 0+2*Y' 4i-2 0)/3, где
Y 4i-3 0 решение в i-3 точке,
Y' 4i 0,Y' 4i-1 0,Y' 4i-2 0 - значения производных в точках i, i-1,
i-2 соответственно.
Для улучшения решения используется корректирующая формула
C 4i+1 0=[9*Y 4i 0-Y 4i-2 0+3*h*(M' 4i+1 0+2*Y' 4i 0-Y' 4i-1 0)]/8, где
M 4i+1 0=P 4i+1 0-112*(P-C 4i 0)/121; M' 4i+1 0=F(X 4i+1 0,M 4i+1 0).
Решение системы в i+1 точке находится по формуле
G=W 4j 0*|P 4i+1,j 0-C 4i+1 0|, где
W 4j 0=1
j- номер компоненты вектора.
На участке "разгона" значения Y 4i-k 0 и Y' 4i-k 0 (k=0, 1, 2)
вычисляются методом Рунге-Кутта по формуле
Y 4i 0=U 4i 5(2) 0-(U 4i 5(i) 0-U 4i 5(2) 0)/15, где i- номер точки, в которой
ищется решение, U 4i 0- решение системы в i-ой точке, полученное с
шагом h/l;
U 5(l) 4i-m+1/l 0=A 5(l) 4i-m 0+(K 5(l) 41 0+2*K 5(l) 42 0+2*K 5(l) 43 0+K 5(l) 44 0) 4i--m+1/l 0/6 4,
где
j=1, 2, ..., n,
K 5(l) 41 0=h*F(X 4i-m 0,A 5(l) 4i-m 0)/l;
K 5(l) 42 0=h*F(X 4i-m 0+h/2*l,A 5(l) 4i-m 0+K 5(l) 41 0/2)/l;
K 5(l) 43 0=h*F(X 4i-m 0+h/2*l,A 5(l) 4i-m 0+K 5(l) 42 0/2)/l;
K 5(l) 44 0=h*F(X 4i-m 0+h/l,A 5(l) 4i-m 0+K 5(l) 43 0)/l.
A, U ,K - векторы n-го порядка; l=1, 2; m=1 при l=1; m=1,
1/2 при l=2;
A 5(l) 4i-1 0=Y 5(l) 4i-1 0; A 5(2) 4i-1/2 0=U 5(2) 4i-1/2 0.
2Характеристика программы. 0
Программа состоит из стандартной информативы, реализующей
описанный метод, рабочей информативы, задающей правые части
уравнений системы и директивы.
Длина стандартной информативы 1600 символов. Объем исход-
ных данных : 7 чисел, 2 массива, n функций. В результате рабо-
ты программы на печать выводится на участке "разгона" X, зна-
чения функций и производных, далее X, G и Y[n] на всем отрезке
интегрирования через Ю шагов и в конце отрезка.
Программа рекомендуется для решения систем обыкновенных
дифференциальных уравнений на больших отрезках, так как счита-
ет быстрее одноточечных методов. Для контроля постоянно выво-
дится погрешность вычислений G, которая позволяет следить за
точностью решения.
"Разгон" (нахождение значений функций и производных в
точках X 40 0, X 40 0+Q, X 40 0+2*Q , X 40 0+3*Q, где Q - шаг интегрирования )
осуществляется методом Рунге-Кутта с увеличенной разрядностью.
В программе предусмотрена возможность при получении боль-
шой погрешности вычисления в точка "разгона" уменьшить шаг ин-
тегрирования в этих точках (см. способ задания J), а при быст-
ром возрастании погрешности вычислений G уменьшить шаг интег-
рирования методом Хемминга или увеличить разрядность вычисле-
ний.
Программа позволяет производить интегрирование как с по-
ложительным, так и с отрицательным шагом (соответственно меня-
ются X 40 0, X 4к 0 и Q).
2Порядок решения задачи. 0
Для решения задачи вводятся стандартная и рабочая инфор-
мативы и директива.
В рабочей информативе после метки Ц программа вычисления
правых частей системы. Здесь 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.
2Описание работы программы
Данная расчетно-графическая работа (далее РГР) составлена
на языке 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. 1 _ 0Hem.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 - векторы точек, для которых производится вычисление.
Для удобства отладки и описания, программа разбита на
части, обозначенные русскими заглавными буквами (Ш,Щ,Л и
т.д.), которые соответствуют блокам, обозначенным в примере
программы, приведенной в задании. Несмотря на то, что приве-
денная программа написана на условном языке, прокомментировать
текст нашей программы на языке MathLab довольно удобно с ис-
пользованием данных обозначений (Конечно, часть блоков опуще-
на, в связи с отсутствием принципиальной значимости. Кроме то-
го изменен порядок появления блоков в программе):
"Э" - начальное вычисление левых частей.
"Ф" - общий цикл, в котором и происходят собственно все
вычисления. Он начинается с конструкций:
While (xk-x1)>0
While 1,
то есть пока не будет достигнут конец, все вычисления
происходят внутри этого цикла.
Также внутри блока "Ф" происходят вычисления корректирую-
щей формулы IAR(i) а также оценка погрешности вычислений G,
если они еще не были рассчитаны на предыдущем шаге.
"Ц" - вычисление левых частей.
"Щ" - на этом этапе происходит перемещение данных в рабо-
чих массивах и X=X+H, то есть происходит переход к следующему
шагу. Также на этом этапе происходит вывод на экран и формиро-
вание выходных массивов Yout, Xout, DGout.
"Л" - в этом блоке происходит расчет самой предсказываю-
щей формулы метода Хемминга - P(i) и Y(i) и происходит расчет
левых частей, т.е. снова в программе появляется блок "Ц".
Затем опять продолжается блок "Ф". Идет проверка на каком
шаге мы находимся и, если он (шаг) меньше 5, то идет подготов-
ка к расчету методом Р-К, то есть идет участок "разгона". Со-
ответственно идет расчет коэффициентов K1, K2, K3 и K4, необ-
ходимых для метода Р-К. Также внутри данного блока еще раз
встречается блок "Ц", в котором происходит расчет левых час-
тей, необходимых для метода Р-К.
Далее происходит проверка шага и уменьшение или увеличе-
ние его в соответствии с заданной точностью. Для возможности
"отката" в случае большого или маленького шага используется
переменная Х1. Также еще раз встречается блок "Ц". Затем, в
случае если все коэффициенты К1-К4 вычислены и шаг удовлетво-
ряет требованиям точности, то происходит расчет методом Р-К,
а также, естественно происходит формирование выходных массивов
Yout, Xout и DGout, а также происходит переход к следующему
шагу (то есть X=X+H) и переход к блоку "Э".
На этом кончается блок "Ф" и вся функция. В начале блока
"Ф" происходит проверка на конец вычислений и если расчеты за-
кончились, то есть мы достигли Xk то происходит возврат в го-
ловную программу. Все выходные данные формируются внутри блока
"Ф", поэтому никаких дополнительных действий не производится.
2Сравнительный анализ и оценка быстродействия
Для сравнения полученных результатов с другими методами
используется метод Адамса, разработанный другой бригадой.
Число операций в методе Хемминга: порядка 2200.
Быстродействие: порядка 0,8 секунды.
Число операций в методе Адамса: порядка 560.
Быстродействие: порядка 0,55 секунды.
(Вычисления проводились на компьютере i486DX4-100)
Как видно из вышеприведенных данных, метод Хемминга проигрывает
по временным показателям и по затратам машинного времени методу Адам-
са, однако стоит заглянуть в приложение, где приведены распечатки гра-
фиков решений и ошибок обоих методов и сразу видно, что метод Адамса
не справляется с контрольным примером для нашей системы, так как ошиб-
ка у него к концу вычислений (Xk=1) возрастает, а в "нашем" методе -
стремится к 0.
2Выводы
Данная РГР по предмету "Численные методы в экономике" ре-
ализует метод Хемминга, который предназначен для решения зада-
чи Коши для обыкновенных дифференциальных уравнений. Програм-
ма, написанная на языке MathLab, хотя и не является оптималь-
ной, но решает поставленную задачу и решает ОДУ довольно боль-
ших степеней сложности, с которыми не справляются другие мето-
ды (например метод Адамса). Это связано с тем, что метод Хем-
минга является многоточечным, а в связи с этим и повышается
точность вычислений, а также устойчивость метода. Однако дан-
ный метод требует больших вычислительных затрат, что связано с
довольно громоздкими формулами а также с большим объемом вы-
числений и поэтому для относительно простых систем целесооб-
разно использовать более простые методы решения.
2Список литературы
1. Д.Мак-Кракен, У.Дорн. "Численные методы и программиро-
вание на Фортране", Издательство "Мир", М. 1977г.
2. О.М.Сарычева. "Численные методы в экономике. Конспект
лекций", Новосибирский государственный технический универси-
тет, Новосибирск 1995г.
3. Н.С.Бахвалов. "Численные методы", Издательство "Нау-
ка", М. 1975г.