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

     Программа предназначена для  численного  решения  системы
обыкновенных дифференциальных уравнений вида:
     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г.