Программа для численного решения системы обыкновенных дифференциальных уравнений 

 

Теоретическая часть
Программа предназначена для численного решения системы обыкновенных дифференциальных уравнений вида:
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 - векторы точек, для которых производится вычисление.
Для удобства отладки и описания, программа разбита на части, обозначенные русскими заглавными буквами (Ш,Щ,Л и т.д.), которые соответствуют блокам, обозначенным в примере программы, приведенной в задании. Несмотря на то, что приведенная программа написана на условном языке, прокомментировать текст нашей программы на языке 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 то происходит возврат в головную программу. Все выходные данные формируются внутри блока "Ф", поэтому никаких дополнительных действий не производится.
Сравнительный анализ и оценка быстродействия
Для сравнения полученных результатов с другими методами используется метод Адамса, разработанный другой бригадой.
Число операций в методе Хемминга: порядка 2200.
Быстродействие: порядка 0,8 секунды.
Число операций в методе Адамса: порядка 560.
Быстродействие: порядка 0,55 секунды.
(Вычисления проводились на компьютере i486DX4-100)
Как видно из вышеприведенных данных, метод Хемминга проигрывает по временным показателям и по затратам машинного времени методу Адамса, однако стоит заглянуть в приложение, где приведены распечатки графиков решений и ошибок обоих методов и сразу видно, что метод Адамса не справляется с контрольным примером для нашей системы, так как ошибка у него к концу вычислений (Xk=1) возрастает, а в "нашем" методе -стремится к 0.
Выводы
Данная РГР по предмету "Численные методы в экономике" реализует метод Хемминга, который предназначен для решения задачи Коши для обыкновенных дифференциальных уравнений. Программа, написанная на языке MathLab, хотя и не является оптимальной, но решает поставленную задачу и решает ОДУ довольно больших степеней сложности, с которыми не справляются другие методы (например метод Адамса). Это связано с тем, что метод Хемминга является многоточечным, а в связи с этим и повышается точность вычислений, а также устойчивость метода. Однако данный метод требует больших вычислительных затрат, что связано с довольно громоздкими формулами а также с большим объемом вычислений и поэтому для относительно простых систем целесообразно использовать более простые методы решения.
Список литературы
1. Д.Мак-Кракен, У.Дорн. "Численные методы и программирование на Фортране", Издательство "Мир", М. 1977г.
2. О.М.Сарычева. "Численные методы в экономике. Конспект лекций", Новосибирский государственный технический университет, Новосибирск 1995г.
3. Н.С.Бахвалов. "Численные методы", Издательство "Наука", М. 1975г.