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

        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)

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