Задания по численным методам

Теоретическая часть

Программа предназначена для численного решения системы

обыкновенных дифференциальных уравнений вида:

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

...

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г.