Неявные методы решения системы уравнений ОДУ

Неявные методы решения системы ОДУ 
 
 

Курсовая  работа 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

ОГЛАВЛЕНИЕ 

ВВЕДЕНИЕ

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

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

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

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

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

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

ЗАКЛЮЧЕНИЕ

СПИСОК  ЛИТЕРАТУРЫ        
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

ВВЕДЕНИЕ 

      Системы дифференциальных уравнений,  зависимости от своей структуры могут быть решены различными методами. Точное решение  получить очень часто не удается,  поэтому мы рассмотрим численные методы решения таких систем.  Далее мы представим две группы  методов: явные и неявные. Для решения систем ОДУ неявными методами придется решать системы нелинейных уравнений,  поэтому придется ввести в рассмотрение  группу  методов решения систем нелинейных уравнений,  которые в свою очередь будут представлены двумя методами. Далее следуют теоретические аспекты описанных методов,  а затем будут представлены описания программ. Сами программы, а также их графики приведены в приложении.

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

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

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

      
 
 
 
 
 

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

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

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

         Неявный метод Эйлера используется  для интегрирования  " жестких  " систем. "Жесткие" системы  это такие системы, в которых 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.  НЕЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА 

         Метод Эйлера  является  методом   Рунге-Кутта  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).