Численные методы в экономике

  

 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. 

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