"Комплект" заданий по численным методам

                                2ВВЕДЕНИЕ

     В экономике очень часто используется модель,  называемая  "черный

ящик",  то есть система у которой известны входы и выходы,  а то,  что

происходит внутри - неизвестно. Законы по которым происходят изменения

выходных  сигналов  в зависимости от входных могут быть различными,  в

том числе это могут быть и дифференциальные законы. Тогда встает проб-

лема решения систем дифференциальных уравнений,  которые в зависимости

от своей структуры могут быть решены различными методами. Точное реше-

ние  получить очень часто не удается,  поэтому мы рассмотрим численные

методы решения таких систем.  Далее мы представим две группы  методов:

явные и неявные. Для решения систем ОДУ неявными методами придется ре-

шать системы нелинейных уравнений,  поэтому придется ввести в рассмот-

рение  группу  методов решения систем нелинейных уравнений,  которые в

свою очередь будут представлены двумя методами. Далее следуют теорети-

ческие аспекты описанных методов,  а затем будут представлены описания

программ. Сами программы, а также их графики приведены в приложении.

     Также стоит отметить, что в принципе все численные методы так или

иначе сводятся к матричной алгебре,  а в экономических  задачах  очень

часто матрицы  имеют  слабую заполненность и большие размеры и поэтому

неэффективно работать с полными матрицами.  Одна из технологий, позво-

ляющая разрешить  данную  проблему  -  технология  разреженных матриц.

В связи с этим, мы рассмотрим данную технологию и операции умножения и

транспонирования над такими матрицами.

     Таким образом мы рассмотрим весь спектр лабораторных работ.  Опи-

сания всех  программ приводятся после теоретической части.  Все тексты

программ и распечатки графиков приведены в приложении.

                          2ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

              1. ОПИСАНИЕ МЕТОДОВ ИНТЕГРИРОВАНИЯ СИСТЕМ ОДУ

                           И ИХ ХАРАКТЕРИСТИК

                  ЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ

        Алгоритм этого метода определяется формулой:

                    x 5m+1 0 = x 5m 0 + h 4m 0*F(x 5m 0, t 4m 0) 4, 0                    (1)

  которая получается путём аппроксимации ряда Тейлора до членов перво-

  го порядка производной x'(t 4m 0),т.к. порядок точности равен 1 (s=1).

        Для аналитического исследования свойств метода Эйлера линеари-

  зуется исходная система ОДУ  x' = F(x, t)  в точке (x 5m 0,t 4m 0):

                         x' = A*x,                               (2)

  где матрица А зависит от точки линеаризации (x 5m 0,t 4m 0).

        Входной сигнал при линеаризации является неизвестной  функцией

  времени и  при  фиксированном t 4m 0 на шаге h 4m 0 может считаться констан-

  той. Ввиду того ,что для линейной системы свойство устойчивости  за-

  висит лишь от А, то входной сигнал в системе (2) не показан. Элемен-

  ты матрицы А меняются с изменением точки линеаризации,т.е. с измене-

  нием m.

        Характеристики метода:

        1.  _Численная устойчивость ..

        Приведя матрицу А к диагональной форме : А = Р* 7l 0*Р 5-1 0 и перейдя

  к новым переменным   y = P 5-1 0*x , система (2) примет вид :

                            y' =  7l 0*y                               (3)

        Тогда метод Эйлера для уравнения (3) будет иметь вид :

                         y 5m+1 0 = y 5m 0 + h* 7l 0 * y 5m 0,                     (4)

  где величина h* 7l 0 изменяется от шага к шагу.

        В этом случае характеристическое уравнение для дискретной сис-

  темы (4) будет иметь вид :

                        1 + h* 7l 0 - r = 0.

  А корень характеристического уравнения равен:

                        r = 1+ h* 7l 0,

  где  7l 0 = 7 a 0  _+ . 7 b 0 .

         _Случай 1 .. Исходная система (3) является асимптотически устой-

  чивой , т.е. нулевое состояние равновесия системы асимптотически ус-

  тойчиво и  7 a 0 < 0.

        Областью абсолютной устойчивости метода интегрирования системы

  (4) будет круг радиусом ,  равным 1 ,  и с центром в точке (0,  -1).

  Таким образом , шаг h должен на каждом интервале интегрирования под-

  бираться таким образом,  чтобы при этом не покидать область А  .  Но

  в таком случае должно выполняться требование :

                            h < 2* 7t 0,                               (5)

  где  7t 0=1/ 72a2 0  - постоянная времени системы (3) .  Она определяет ско-

  рость затухания переходных процессов в ней. А время переходного про-

  цесса T 4пп 0 = 3* 7t 0 , где  7t 0 =  72a2 5-1 0.

        Если иметь ввиду, что система (3) n-го порядка, то

                         T 4пп 0 > 3* 7t 4max 0,

  где  7t 4max 0 =  72a 4min 72 5-1 7  0;  7a 4min  0= min  7a 4i 0 . Из всех  7a 4i 0 (i=[1;n]) выбирает-

  ся максимальное значение : max 72a 4i 72 0= 7a 4max 0  и тогда можно вычислить :

                         7t 4min  0= 1/ 7a 4max 0,

  а шаг должен выбираться следующим образом :

                   h < 2/ 7a 4max 0  или   h < 2* 7t 4min 0.

         _Случай 2 ..  Нулевое состояние равновесия системы (2) неустойчи-

  во, т.е.  7a 0 > 0 . Т.е. система тоже неустойчива , т.е.  72 0r 72 0>1. Поэтому

  областью относительной  устойчивости  явного  метода Эйлера является

  вся правая полуплоскость . Выполняется требование :

                     72 0 1+h* 7l 0  72  0< 7 2  0e 5hl 7 2 0                            (6)

        2.  _Точность метода ..

        Так как  формула интегрирования (1) аппроксимирует ряд Тейлора

  для функции x(t 4m 0+1) до линейного по h члена включительно. Существует

  такое значение t в интервале [t 4m 0 , t 4m+1 0], при котором

                      Е 4i 5am 0 =1/2! * h 4m 52 0*x 4i 0''(t),

  где i=[1;n].

        3.  _Выбор шага интегрирования ..

        Должны соблюдаться условия абсолютной  (5)  или  относительной

  (6) устойчивости в зависимости от характера интегрируемой системы.

        Для уравнения первого порядка шаг должен быть :

                              h < 2* 7t 0 .

        Для уравнений n-ого порядка :

                             h 4i 0 <= 2* 7t 4i  0.

  А окончательно шаг выбирают по условиям устойчивости :

                    h < 2* 7t 4min 0 ,   7t 4min 0 = min  7t 4i

  Вначале задаётся допустимая ошибка аппроксимации ,  а в процессе ин-

  тегрирования шаг подбирается следующим образом :

        1) по формуле (1) определяется очередное значение x 5m+1 0;

        2) определяется dx 4i 5m 0 = x 4i 5m+1 0 - x 4i 5m 0 ;

        3) условие соблюдения точности имеет вид :

                h 4i 5m 0 <= E 4i 5aдоп 7/ 0[ 72 0f 4i 0(x 5m 0,t 4m 0) 72 0 + E 4i 5aдоп 7/ 0h 4max 0]           (7)

        4) окончательно на m-м интервале времени выбирается в виде:

                              h 4m 0 = min h 4i 5m 0.

                       ЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА

        Метод Эйлера  является  методом  Рунге-Кутта  1-го  порядка  .

  Методы Рунге-Кутта  2-го  и  4-го  порядка  являются  одношаговыми ,

  согласуются с рядом Тейлора до порядка точности s  ,  который  равен

  порядку метода  .  Эти  методы  не  требуют  вычисления  производных

  функций , а только самой функции в нескольких точках на шаге h 4m 0.

        Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем :

                x 5m+1 0 = x 5m 0 + h 4m 0/2 (k 41 0+k 42 0),

  где k 41 0=f(x 5m 0,t 4m 0) ; k 42 0=f(x 5m 0+h 4m 0*k 41 0,t 4m 0+h 4m 0).

        Ошибка аппроксимации Е 5a 0 = k*h 4m 53 0 .

        Алгоритм метода Рунге-Кутта 4-го порядка

                x 5m+1 0=x 5m 0+h 4m 0/6(k 41 0+2k 42 0+2k 43 0+k 44 0),

  где k 41 0=f(x 5m 0,t 4m 0); k 42 0=f(x 5m 0+h 4m 0/2*k 41 0,t 4m 0+h 4m 0/2); k 43 0=f(x 5m 0+h 4m 0/2*k 42 0,t 4m 0+h 4m 0/2);

                       k 44 0=f(x 5m 0+h 4m 0*k 43 0,t 4m 0+h 4m 0).

        Ошибка аппроксимации Е 5a 0 = k*h 4m 55 0.

              НЕЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ

         Неявный метод Эйлера используется для интегрирования  " жест-

  ких " систем. "Жесткие" системы это такие системы, в которых 7  0 ( 7l 4max 0)

  и ( 7l 4min 0) сильно отключаются друг от друга , то в решениях системы

                           x' = A*x                              (1)

  будут присутствовать экспоненты,  сильно отличаются друг от друга по

  скорости затухания .  Шаг интегрирования для таких систем должен вы-

  бираться по условиям устойчивости из неравенства

                          h <= 2* 7t 4min , 0                          (2)

  где  7t 0=1/ 72a2 0  - постоянная времени системы  y' =  7l 0*y . Она определяет

  скорость затухания  переходных  процессов  в  ней .  Неравенство (2)

  должно выполняться на всем участке решения , что соответственно тре-

  бует значительных затрат машинного времени.

         Алгоритм этого метода определяется формулой:

                    x 5m+1 0 = x 5m 0 + h 4m 0*F(x 5m+1 0, t 4m+1 0)  4                 0(3)

         Если h 4m 0 мал, то x 5m 0 и х 5m+1 0 близки друг к другу. В качестве на-

  чального приближения берется точка x 5m 0 , а следовательно , между x 5m 0 и

  x 5m+1 0 будет существовать итерационный процесс.

        Для аналитического  исследования свойств  метода Эйлера линеа-

  ризуется исходная система ОДУ  x' = F(x, t)  в точке (x 5m 0,t 4m 0):

                         x' = A*x,

  где матрица А зависит от точки линеаризации (x 5m 0,t 4m 0).

        Входной сигнал при линеаризации является неизвестной  функцией

  времени и  при  фиксированном t 4m 0 на шаге h 4m 0 может считаться констан-

  той. Ввиду того ,что для линейной системы свойство устойчивости  за-

  висит лишь от А, то входной сигнал в системе (1) не показан. Элемен-

  ты матрицы А меняются с изменением точки линеаризации,т.е. с измене-

  нием m.

        Характеристики метода:

        1.  _Численная устойчивость ..

        Приведя матрицу А к диагональной форме : А = Р* 7l 0*Р 5-1 0 и перейдя

  к новым переменным   y = P 5-1 0*x , система (3) примет вид :

                            y' =  7l 0*y                             (4)

        Тогда метод Эйлера для уравнения (4) будет иметь вид :

                         y 5m+1 0 = y 5m 0 + h* 7l 0 * y 5m+1 0,                 (5)

  где величина h* 7l 0 изменяется от шага к шагу.

        В этом случае характеристическое уравнение для дискретной сис-

  темы (5) будет иметь вид :

                        1 - h* 7l 0*r - 1 = 0.

  А корень характеристического уравнения равен:

                         r = 1/(1-h* 7l 0) ,

  где  7l 0 = 7 a 0  _+ . 7 b 0 .

         _Случай 1 .. Исходная система (4) является асимптотически устой-

  чивой , т.е. нулевое состояние равновесия системы асимптотически ус-

  тойчиво и  7 a 0 < 0.

        Областью абсолютной устойчивости метода интегрирования системы

  (5) будет вся левая полуплоскость. Таким образом , шаг  h должен  на

  каждом интервале интегрирования подбираться таким образом, чтобы при

  этом не покидать эту область.  Но в таком случае должно  выполняться

  требование :

                            h < 2* 7t 0,                             (6)

  где  7t 0=1/ 72a2 0  - постоянная времени системы (4) .  Она определяет ско-

  рость затухания переходных процессов в ней. А время переходного про-

  цесса T 4пп 0 = 3* 7t 0 , где  7t 0 =  72a2 5-1 0.

        Если иметь ввиду, что система (4) n-го порядка, то

                         T 4пп 0 > 3* 7t 4max 0,

  где  7t 4max 0 =  72a 4min 72 5-1 7  0;  7a 4min  0= min  7a 4i 0 . Из всех  7a 4i 0 (i=[1;n]) выбирает-

  ся максимальное значение : max 72a 4i 72 0= 7a 4max 0  и тогда можно вычислить :

                         7t 4min  0= 1/ 7a 4max 0,

  а шаг должен выбираться следующим образом :

                   h < 2/ 7a 4max 0  или   h < 2* 7t 4min 0.

         _Случай 2 ..  Нулевое состояние равновесия системы (4) неустойчи-

  во, т.е.   7a 0  >  0  .  Т.е.  система тоже неустойчива ,  т.е.   72 0r 72 0>1,

  а следовательно :

                          72 0 1/(1-h* 7l 0)  72 0 > 1.

         С учетом ограничения на скорость изменения приближенного  ре-

  шения относительно точного :

                      72 0 1/(1-h* 7l 0)  72 0 > 7 2  0e 5hl 7 2 0.                    (7)

         Из этого соотношения следует , что при  72 0h* 7l2 0 -> 1 левая часть

  стремится к бесконечности . Это означает , что в правой полуплоскос-

  ти есть некоторый круг радиусом , равным 1 , и  с  центром  в  точке

  (0, 1), где неравенство (7) не выполняется.

         2.  _Точность метода ..

         Ошибка аппроксимации  по  величине равна ошибке аппроксимации

  явного метода Эйлера , но она противоположна по знаку :

                      Е 4i 5am 0 =-1/2! * h 4m 52 0*x 4i 0''(t),

  где t 4m 0 <= t <= t 4m+1 0,

      i=[1;n].

         3.  _Выбор шага интегрирования ..

         Должны соблюдаться условия абсолютной (6)  или  относительной

  (7) устойчивости в зависимости от характера интегрируемой системы.

         Для уравнения первого порядка шаг должен быть :

                              h < 2* 7t 0 .

         Для уравнений n-ого порядка :

                             h 4i 0 <= 2* 7t 4i  0.

         Однако область абсолютной устойчивости - вся левая  полуплос-

  кость . Поэтому шаг с этой точки зрения может быть любым.

         Условия относительной устойчивости жестче,  так как есть  об-

  ласть , где они могут быть нарушены. Эти условия будут выполняться ,

  если в процессе решения задачи контролировать ошибку аппроксимации ,

  а шаг корректировать .  Таким образом, шаг можно выбирать из условий

  точности, при этом условия устойчивости будут соблюдены автоматичес-

  ки. Сначала задается допустимая погрешность аппроксимации :

                    E 4i 5aдоп 0 <= 0,001  72 0x 4i 72 4max 0,

  где i=[1;n].

         Шаг выбирается в процессе интегрирования следующим образом:

         1. Решая систему (3) относительно x 5m+1 0 с шагом  h 4m 0,  получаем

  очередную точку решения системы x = F(x,t) ;

         2. Оцениваем величину ошибки аппроксимации по формуле:

      Е 4i 5am 0 =   72 0h 4m 7/ 0(h 4m 0+h 4m-1 0)*[(x 4i 5m+1 0  - x 4i 5m 0) - h 4m 7/ 0h 4m-1 0*(x 4i 5m 0 -x 4i 5m-1 0)] 72

         3. Действительное значение аппроксимации сравнивается  с  до-

  пустимым. Если Е 4i 5am 0 < E 4i 5aдоп 0, то выполняется следующий пункт, в про-

  тивном случае шаг уменьшается в два раза ,  и вычисления повторяются

  с п.1.

         4. Рассчитываем следующий шаг:

                h 4i 5m+1 0 = SQR(E 4i 5aдоп 7/2 0Е 4i 5am 72 0) * h 4m 0.

         5. Шаг выбирается одинаковым для всех элементов вектора X :

                              h 4m+1 0 = min h 4i 5m+1 0.

         6. Вычисляется новый момент времени и алгоритм повторяется .

                       НЕЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА

         Метод Эйлера  является  методом  Рунге-Кутта  1-го  порядка .

  Методы Рунге-Кутта  2-го  и  4-го  порядка  являются  одношаговыми ,

  согласуются с рядом Тейлора до порядка точности s  ,  который  равен

  порядку метода  .  Эти  методы  не  требуют  вычисления  производных

  функций , а только самой функции в нескольких точках на шаге h 4m 0.

         Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем:

                x 5m+1 0 = x 5m 0 + h 4m 0/2 (k 41 5m+1 0+k 42 5m+1 0),

  где    k 41 5m+1 0=f(x 5m+1 0,t 4m+1 0);  k 42 5m+1 0=f(x 5m+1 0-h 4m 0*k 41 5m+1 0,t 4m+1 0).

         Ошибка аппроксимации Е 4m 5a 0 = k*h 4m 53 0 .

         Алгоритм метода Рунге-Кутта 4-го порядка

                x 5m+1 0 = x 5m 0 + h 4m 0/6 (k 41 5m+1 0+2k 42 5m+1 0+2k 43 5m+1 0+k 44 5m+1 0),

  где    k 41 0=f(x 5m+1 0,t 4m+1 0);     k 42 0=f(x 5m+1 0-h 4m 0/2*k 41 5m+1 0,t 4m+1 0-h 4m 0/2);

         k 43 0=f(x 5m+1 0-h 4m 0/2*k 42 5m+1 0,t 4m+1 0-h 4m 0/2);     k 44 0=f(x 5m+1 0-h 4m 0*k 43 5m+1 0,t 4m 0).

         Ошибка аппроксимации Е 5a 0 = k*h 4m 55 0.

                   2. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ САУ

                            МЕТОД НЬЮТОНА

     Итерационная формула метода Ньютона имеет вид

                     X 5m+1 0=X 5m  0- 5  0J 5-1  0* 5  0(X 5m 0) 5  0* 5  0F(X 5m 0),

     где J(X)=F 4x 5| 0/ 4x=xm

     Характеристики метода:

     1. Сходимость. Покажем, что в точке P(G 4x 5| 0(X 5* 0))=0

     Здесь G(x)=X  - J 5-1 0(x) * F(x) - изображение итерационного процес-

са. Условие сходимости метода:  P(G 4x 5| 0(X)) < 1 должно  выполняться  для

всех значений  X 5m 0.  Это  условие осуществляется при достаточно жестких

требованиях к F(x):  функция должна быть непрерывна и  дифференцируема

по X, а также выпукла или вогнута вблизи X 5* 0. При этом выполняется лишь

условие локальной сходимости.  Причем можно показать,  что чем ближе к

X 5* 0, тем быстрее сходится метод.

     2. Выбор начального приближения X 50 0 - выбирается достаточно близко

к точному  решению.  Как именно близко - зависит от скорости изменения

функции F(x) вблизи X 5* 0:  чем выше скорость,  тем меньше  область,  где

соблюдается условие  сходимости и следовательно необходимо ближе выби-

рать X 50 0 к X 5* 0.

     3. Скорость сходимости определяется соотношением

                  ║E 5m+1 0║ = C 4m 0 * ║E 5m 0║ 51+p 0,  0 < P < 1.

     При X -> X 5* 0 величина P -> 1. Это значит, что вблизи точного реше-

ния скорость сходимости близка к квадратичной.  Известно,  что  метода

Ньютона сходится на 6-7 итерации.

     4. Критерий окончания итераций - здесь могут использоваться  кри-

терии точности, то есть максимальная из норм изменений X и F(x).

                       ДИСКРЕТНЫЙ МЕТОД НЬЮТОНА

     Трудность использования метода Ньютона состоит в

     1. Необходимости вычисления на каждом этапе J=F 4x 5| 0.

     Возможно несколько путей решения:

     - аналитический способ. Наиболее эффективен, однако точные форму-

лы могут быть слишком большими, а также функции могут быть заданы таб-

лично.

     - конечно-разностная аппроксимация:

       dF 4i 0   F 4i 0(x 41 0,...,x 4j 0+dx 4j 0,...,x 4n 0) - F 4i 0(x 41 0,...,x 4j 0-dx 4j 0,...x 4n 0)

       ─── = ──────────────────────────────────────────────────

       dX 4j 0                           2*dX 4j

     В этом случае мы имеем уже дискретный метод Ньютона,  который уже

не обладает квадратичной сходимостью. Скорость сходимости можно увели-

чить, уменьшая dX 4j 0 по мере приближения к X 5* 0.

     2. Вычисление матрицы J 5-1 0 на каждом шаге требует значительных вы-

числительных затрат, поэтому часто решают такую систему:

                         dX 5  0= 5  0J 5-1 0(X 5m 0) 5  0* 5  0F(X 5m 0)

     или умножая правую и левую часть на J, получим:

                         J(X 5m 0) 5  0* 5  0dX 5m  0= 5  0F(X 5m 0)

     На каждом M-ом шаге матрицы J и F известны. Необходимо найти dX 5m 0,

как решение вышеприведенной системы линейных АУ, тогда

                             X 5m+1 0=X 5m 0+dX 5m

     Основной недостаток  метода  Ньютона  - его локальная сходимость,

поэтому рассмотрим еще один метод.

                МЕТОД ПРОДОЛЖЕНИЯ РЕШЕНИЯ ПО ПАРАМЕТРУ

     Пусть t - параметр,  меняющийся от 0 до 1.  Введем в рассмотрение

некоторую систему

                              H(X,t)=0,

     такую, что:

     1. при t=0 система имеет решение X 50

     2. при t=1 система имеет решение X 5*

     3. вектор-функция H(X,t) непрерывна по t.

     Вектор функция может быть выбрана несколькими способами, например

                        H(X,t) = F(X) + (t-1)

                                 или

                          H(X,t) = t * F(X)

     Нетрудно проверить, что для этих вектор-функций выполняются усло-

вия 1-3.

     Идея метода состоит в следующем:

     Полагаем t 41 0= 7D 0t и решаем систему H(X,t 41 0)=0 при выбранном X 50 0. Полу-

чаем вектор  X 5t1 0.  Далее берем его в качестве начального приближения и

решаем при новом t 42 0=t 41 0+ 7D 0t систему H(X,t 42 0)=0,  получаем X 5t2 0 и так далее

до тех пор, пока не будет достигнута заданная точность.

                   3. ТЕХНОЛОГИЯ РАЗРЕЖЕННЫХ МАТРИЦ

                        ОСНОВНЫЕ ИДЕИ МЕТОДА.

     Основные требования к этим методам можно свести в 3 утверждения:

     1. Хранить в памяти только ненулевые элементы.

     2. В процессе решения принимать меры, уменьшающие возможность по-

явления новых ненулевых элементов.

     3. Вычисления производить только с ненулевыми элементами.

                     Разреженный строчный формат

     Это одна  из  широко  используемых  схем для хранения разреженных

матриц, которая предъявляет минимальные требования к  памяти  и  очень

удобна для выполнения основных операций с матрицами.

     Значения ненулевых элементов и соответствующие столбцовые индексы

хранятся по  строкам  в двух массивах:  NA и JA.  Для разметки строк в

этих массивах используется массив указателей  IA,  отмечающих  позиции

массивов AN и JA, с которых начинается описание очередной строки. Пос-

ледняя цифра в массиве IA содержит указатель первой свободной  позиции

в JA и AN.  Рассмотрим конкретный пример,  который будет также и далее

применятся для демонстрации других операций и который был  использован

в качестве  контрольного  примера для программы,  выполняющей основные

операции с разреженными матрицами.

          ┌           ┐

          │ 3 0 0 5 0 │      Позиция:  1 2  3 4  5 6  7 8  9 10

          │ 0 4 0 0 1 │           IA:  1    3    5    7    9    11

      A = │ 0 0 8 2 0 │           JA:  1 4  2 5  3 4  1 4  2 5

          │ 5 0 0 6 0 │           AN:  3 5  4 1  8 2  5 6  7 9

          │ 0 7 0 0 9 │

          └           ┘

     Данный способ представления называется полным (так как  представ-

лена вся  матрица  А)  и упорядоченным (так как элементы каждой строки

хранятся в соответствии с возрастанием столбцовых индексов). Обознача-

ется RR(c)O - Row-wise representation Complete and Ordered (англ.).

     Массивы IA и JA представляют портрет матрицы А.  Если в алгоритме

разграничены этапы символической и численной обработки,  то массивы IA

и JA заполняются на первом этапе, а массив AN - на втором.

                 Транспонирование разреженной матрицы

     Пусть IA,  JA,  AN - некоторое представление разреженной матрицы.

Нужно получить IAT,  JAT, ANT - разреженное представление транспониро-

ванной матрицы. Алгоритм рассмотрим на нашем примере:

     IA = 1 3 5 7 9 11

     JA = 1 4 2 5 3 4 1 4 2 5

     AN = 3 5 4 1 8 2 5 6 7 9

     Символический этап.

     1. Заводим 5 целых списков по числу столбцов матрицы А.  пронуме-

руем их от 1 до 6.

     2. Обходим 1 строку и заносим 1 в те списки,  номера которых ука-

заны в JA:

                        1: 1

                        2:

                        3:

                        4: 1

                        5:

     3. Обходим вторую строку, вставляя аналогичным образом 2:

                        1: 1

                        2: 2

                        3:

                        4: 1

                        5: 2

     В итоге получим:

                        1: 1 4

                        2: 2 5

                        3: 3

                        4: 1 3 4

                        5: 2 5

     Массив ANT можно получить,  вставляя численное  значение  каждого

ННЭ, хранимое в AN,  в вещественный список сразу после того, как соот-

ветствующее целое внесено в целый список. В итоге получим:

     IAT = 1 3 5 6 9 11

     JAT = 1 4 2 5 3 1 3 4 2 5

     ANT = 3 5 4 7 8 5 2 6 1 9

                  Произведение разреженной матрицы и

                     заполненного вектора-столбца

     Алгоритм рассмотрим на нашем конкретном примере:

     IAT = 1 3 5 7 9 11

     JAT = 1 4 2 5 3 1 3 4 2 5

     ANT = 3 5 4 7 8 5 2 6 1 9

     B = ( -5 4 7 2 6 )

     Обработка 1 строки:

     Просматриваем массив IA и обнаруживаем, что первая строка матрицы

А соответствует  первому  и  второму  элементу  массива  JA:  JA(1)=3,

JA(2)=4, то есть ННЭ являются a 411 0 и a 414 0.

     Просматриваем массив AN и устанавливаем, что a 411 0=3 и a 414 0=5.

     Обращаемся к вектору B,  выбирая из него соответственно  b 41 0=-5  и

b 44 0=2.

     Вычисляем C 41 0= 3 * (-5) + 5 * 2 = -5.

     Абсолютно аналогично, вычисляя остальные строки, получим:

                         C = (-5 58 56 1 58).

                 Произведение двух разреженных матриц

     Рассмотрим метод на конкретном примере.  Предположим, что нам не-

обходимо перемножить две матрицы:

     IA = 1 3 5 7 9 11

     JA = 1 4 2 5 3 4 1 4 2 5

     AN = 3 5 4 1 8 2 5 6 7 9

     IAT = 1 3 5 7 9 11

     JAT = 1 4 2 5 3 1 3 4 2 5

     ANT = 3 5 4 7 8 5 2 6 1 9

     Суть метода состоит в том, что используя метод переменного перек-

лючателя и расширенный вещественный накопитель, изменяется порядок на-

копления сумм для формирования элементов матрицы С.  Мы будем формиро-

вать элементы C 4ij 0 не сразу,  а постепенно накапливая, обращаясь только

к строкам матрицы B.

     Символический этап.

     Определяем мерность IX = 5

     IX = 0 0 0 0 0

     1-я строка матрицы JAT: 1 4

       JA(1) = 1 4      JC(1) = 1 4

       IX = 1 0 0 1 0

       JA(4) = 1 4

       IX(1) = 1 ?  Да. JC(1) - без изменений

       IX(4) = 1 ?  Да. JC(1) - без изменений

     2-я строка матрицы JAT: 2 5

       JA(2) = 2 5      JC(2) = 2 5

       IX = 1 2 0 1 2

       JA(5) = 2 5   -> JC(2) - без изменений

     ....

     4-я строка матрицы JAT: 1 3 4

       JA(1) = 1 4      JC(4) = 1 4

       IX = 4 2 2 4 2

       JA(3) = 3 4

       IX(3) = 4 ? Нет. JC(4) = 1 4 3

       IX(4) = 4 ? Да.  JC(4) - без изменений

     ....

     в итоге получим:

     IC = 1 3 5 7 10 12

     JC = 1 4 2 5 3 4 1 4 3 2 5

     Численный этап.

     X = 0 0 0 0 0

     1-я строка: JC(1) = 1 4

     AN(1) = 3 5,

        AA = 3

           ANT(1) = 3 5

           AA * ANT = 9 15

           X = 9 0 0 15 0

        AA = 5

           ANT(1) = 3 5

           AA * ANT = 15 25

           X = 24 0 0 40 0

     CN(1) = 24 40

     ....

     Аналогично проделывая действия над другими строками получим:

     IC: 1       3       5      7          10     12

     JC: 1  4    2  5    3  4   1  4  3    2  5

     CN: 24 40   20 35   80 0   55 22 66   16 144

     Все вышеприведенные операции были получены на  компьютере  и  ре-

зультаты совпали  для нашего контрольного примера.  Описание программы

на языке 2 C 0,  занимающейся этими операциями не приводится, так как дан-

ная программа нами не разрабатывалась, однако в приложении присутству-

ет распечатка этой программы.

                 2ПРАКТИЧЕСКАЯ ЧАСТЬ. ОПИСАНИЯ ПРОГРАММ.

     1. ЯВНЫЕ МЕТОДЫ.

     Данная группа представлена 3 программами, реализующими методы Эй-

лера,Рунге-Кутта 2-го и 4-го порядков.  В данной группе все  программы

индентичны, поэтому далее следует описание программы,  реализующем ме-

тод Эйлера,  с указанием различий для методов Рунге-Кутта 2-го и  4-го

порядков.

      1EILER.M

     Головной модуль.

     Входные и выходные данные: отсутствуют.

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Стандартный головной модуль.  Происходит очистка экрана,  задание

начальных значений по времени и по Y.  Затем происходит вызов подпрог-

раммы EIL.M (RG2.M или RG4.M для методов Рунге-Кутта 2 и 4 порядков) а

после получения результатов строятся графики.

      1EIL.M

     Вычислительный модуль.

     Входные данные:

     FunFcn -  имя  подпрограммы,  написанной  пользователем,  которая

возвращает левые части уравнения для определенного момента времени.

     T0, Tfinal - начальные и конечные моменты времени.

     Y0 - начальные значения.

     Выходные данные:

     Tout - столбец времени

     Yout - столбцы решений по каждой координате

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Данный модуль  и реализует собственно метод Эйлера (Рунге-Кутта 2

или 4-го порядков). В начале происходит инициализация внутренних пере-

менных, а также вычисление максимального шага, который затем использу-

ется для определения точности.  Далее следует цикл While...End  внутри

которого и выполняются все необходимые действия:  по формуле (для каж-

дого метода своя!) вычисляется значение Y на каждом шаге (при  необхо-

димости вызывается  подфункция  FunFcn)  а  также формируются выходные

значения Tout и Yout. Далее метод сам корректирует свой шаг, по форму-

ле описанной выше (в теоретическом разделе).  Этот цикл выполняется до

тех пор, пока значение Tfinal не будет достигнуто. Все выходные значе-

ния формируются  внутри данного цикла и поэтому никаких дополнительных

действий не требуется. В связи с этим с окончанием цикла заканчивается

и подпрограмма  EIL.M.  Методы  Рунге-Кутта реализованы аналогично,  с

учетом отличий в формулах вычислений (более сложные по сравнению с ме-

тодом Эйлера).

     2. НЕЯВНЫЕ МЕТОДЫ.

     Представлены группой из двух похожих между собой программ, реали-

зующих соответственно неявные методы Эйлера и Рунге-Кутта  2  порядка.

Также как и в вышеприведенном случае, будет описан метод Эйлера, а от-

личия метода Рунге-Кутта будут отмечены в скобках.

      1NME.M

     Головной модуль.

     Входные и выходные данные отсутствуют.

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Выполняет стандартные действия:  очистка экрана,  инициализация и

ввод начальных значений,  вызов подпрограмм обработки и вычислений,  а

также построение графиков.

      1NMEF.M, NRG2.M

     Вычислительные модули.

     Входные данные:

     T0, Tfinal - начальные и конечные моменты времени

     X0 - вектор-столбец начальных значений.

     H - начальный шаг

     A - матрица,  на диагонали которой стоят собственные числа линеа-

ризованной системы ОДУ.

     Выходные данные:

     T - столбец времени

     X - столбец решений

      7D 0X - столбец ошибки

     Пояснения к тексту модуля:

     Стандартные действия:  инициализация  начальных  значений ,  цикл

While T < Tfinal,  вычисление по формулам, вывод промежуточных резуль-

татов, формирование выходных значений массивов.

     3. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ САУ

     Представлены группой из 4-х методов: метод последовательных приб-

лижений, метод Ньютона,  метод Ньютона дискретный,  метод  продолжения

решения по параметру.

                 Метод последовательных приближений.

      1MMPP.M

     Головной модуль.

     Входные и выходные данные отсутствуют.

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Очистка экрана,  инициализация начальных значений, вызов вычисли-

тельного модуля MPP.M, вывод результатов в виде графиков.

      1MPP.M

     Вычислительный модуль.

     Входные данные:

     X0 - начальное приближение в виде вектора-строки

     Fun1 - функция, возвращающая вычисленные левые части

     Fun2 - функция, возвращающая матрицу Якоби в определенной точке.

     E - допустимая ошибка.

     Выходные данные:

     Mout - номера итераций

     Xout - приближения на каждой итерации

     DXout - ошибка на каждой итерации

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Аналогичен вышеприведенным вычислительным модулям - инициализация

начальных значений,  вычисления  по формулам,  вывод промежуточных ре-

зультатов, формирование выходных значений. По мере необходимости вызы-

вает подпрограммы Fun1 и Fun2.

                 Методы Ньютона и Ньютона дискретный

      1MNEWT.M

     Головной модуль.

     Входные и выходные данные отсутствуют.

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Стандартный головной модуль  -  выполняет  действия,  аналогичные

предыдущим головным модулям. Вызывает из себя NEWT.M и NEWTD.M - моду-

ли реализующие методы Ньютона и Ньютона дискретный,  а также строит их

графики на одной координатной сетке.

      1NEWT.M, NEWTD.M

     Вычислительные модули.

     Входные данные:

     X0 - начальное приближение в виде вектора-строки

     Fun1 - функция, возвращающая левые части

     Fun2 - функция,  вычисляющая матрицу Якоби (только для метода

                                                             Ньютона!)

     E - допустимая ошибка

     Выходные данные:

     Mout - номера итераций

     Xout - приближения на каждой итерации

     DXout - ошибка на каждой итерации

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модулей:

     Стандартные вычислительные модули,  производящие  соответствующие

им действия. Отличие их в том, что в первом случае для вычисления мат-

рицы Якоби вызывается подпрограмма,  а во втором случае матрица  Якоби

вычисляется внутри модуля.

                Метод продолжения решения по параметру

      1MMPRPP.M

     Головной модуль.

     Входные и выходные данные отсутствуют.

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Стандартный головной модуль со всеми вытекающими отсюда  последс-

твиями.

      1MPRPP.M

     Вычислительный модуль.

     Входные данные:

     Fun1 - имя подпрограммы, вычисляющей правые части

     Fun2 - имя подпрограммы, вычисляющем матрицу Якоби

     X0 - начальное приближение

     dT - начальный шаг

     Edop - допустимая ошибка

     Trace - вывод на экран

     Язык реализации: PC MathLab

     Операционная система: MS-DOS 3.30 or higher

     Пояснения к тексту модуля:

     Стандартный вычислительный модуль  -  инициализация,  вычисление,

вывод, формирование  результата.  Стоит отметить,  что поскольку метод

имеет глобальную сходимость,  то объем вычислений тут значительный,  а

это значит, что при выполнении вычислений требуется значительное коли-

чество свободной оперативной памяти.

                                 2ВЫВОДЫ

     При выполнении данных лабораторных работ  были  изучены  основные

численные методы для решения ОДУ,  САУ, а также технология разреженных

матриц. Заодно были получены основные навыки в программировании  мате-

матической системы  PC  MathLab.  Каждый  из представленных методов по

своему хорош и применяется для систем определенного вида.