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

                                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.  Каждый  из представленных методов по своему хорош и применяется для систем определенного вида.