• уравнения конечных разностей
  • разностные уравнения в векторной форме
  • forward Euler' method (FEM)
  • backward Euler' method (BEM)
  • метод трапеций
  • метод скачков
  • backward-forward method (BFM)
  • метод Runge-Kutta

  • численное дифференцирование

    уравнения конечных разностей

         t  =  0, 1, 2, ...
    
         s  = ( x₁
                x₂
                .
                xₙ )
    
         s₁ = A * s₀
         s₂ = A * s₁
         ..
    
         sn - sk  =  An-k * sk
    

    разностные уравнения в векторной форме

    старинное уравнение (операция) Ньютона A * x = y , yi = xi+1 - xi

    пусть есть последовательность временных квантов:   t0,  ∆t + t0,  ....,  k * ∆t + t0

      пусть есть вектор x = ( x₀
                              x₁
                              .
                              xₙ )
    
    пусть yt = A * xt

    тогда, полагая ∆t = 1, имеем yt+1 - yt = ∆yt

    и матрица такого преобразования (полагая квант ∆t равным 1):

    
        (  y'₁        ( -1  1                         ( x₀ 
           y'₂             -1  1                        x₁ 
           y'₃    =           -1  1          *  A  *    x₂
           .                     .                      . 
           y'ₙ )                    -1  1 )             xₙ )
    
    

    для разности второго порядка: 2 * yt - yt+1 - yt-1 = ∆∆yt

    и матрица такого преобразования:

    
        ( y"₁       (  2 -1                           ( x₀ 
          y"₂         -1  2 -1                          x₁ 
          y"₃    =       -1  2 -1            *  A  *    x₂ 
          .                  .                          . 
          y"ₙ )                      -1  2 )            xₙ )
    
    

    forward Euler' method (FEM)

         sn+1 - sn  =  A * s n 
    

    матрица преобразования:

         sn+1 = (I  +  A) * sn 
    

    в коде октава это будет что-то вроде того:

    % forward Euler' method 
    %  u" + u = 0;  u(0) = 1,  u'(0) = 0  
    
    %  u' =   v
    %  v' = - u
    %  un+1 = un  +   ∆ * vn
    %  vn+1 = vn  -   ∆ * un
    
    %   |un+1|      | 1  ∆|     |un|
    %   |    |  =   |     |  *  |  |
    %   |vn+1|      |-∆  1|     |vn|
    
    
    function fem (n)
    
      x = y = zeros (n, 1) ;
      d = 2 * pi / 100 ;                 % step
    
      x(1) = 1 ; y(1) = 0 ;              % initial conditions
    
      A = [1 d ; -d 1] ; 
    
      for t = 2 : n
        Z = A * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ;
      end
    
      plot (x, y, 'linewidth', 2)
    
    end
    

    метод быстрый, но нестабильный, что легко видеть на графике - спираль раскручивается. сумма собственных значений равна трейсу матрицы (т.е. 2), а произведение равно детерминанту (т.е. 1 + d². матрица несимметричная, значит собственные значения - комплексные и одно из них по модулю больше 1, откуда и следует нестабильность системы. таково се ля ви

    график на фазовой плоскости

    разница между соседними витками спирали характеризует точность

    ошибка на одном шаге алгоритма составляет u" * (∆t)²/2

    суммарная ошибка за N шагов будет (1 + A * (∆t)²)N * err₁ + .. + (1 + A * (∆t)²)N-k * errk + .. + errn

    откуда следует критерий стабильности решения: |1 + A * (∆t)²| ≤ 1

    в общем случае (не для данного уравнения) если решение стабильно, то суммарная ошибка будет ограничена сверху значением: u"max * N * (∆t)²/2

    заметив, что N * ∆t = T , где T - общее время наблюдения за системой (в тиках), получаем: err ~ T/2 * ∆t * u"max

    backward Euler' method (BEM)

        sn+1 - sn = A *  sn+1
    

    ЗАМЕЧАНИЕ: метод неявный, потому что искомое есть как в lhs, так и в rhs, в отличие от прямого метода, где к моменту вычисления sn+1 значение sn было уже известно (либо из начальных условий на первом шаге, либо как результат предыдущего шага)

    матрица преобразования:

         sn+1 = (I  -  A)⁻ * sn
    

    в октаве можно наваять что-либо в духе:

    % backward Euler' method 
    %  u" + u = 0     u(0) = 1     u'(0) = 0 
    
    %  u' =   v
    %  v' = - u
    %  un+1 - un =   ∆ * vn+1
    %  vn+1 - vn = - ∆ * un+1
    
    %   |un+1 - un|     | 0  ∆|     |un+1|
    %   |         |  =  |     |  *  |    |
    %   |vn+1 - vn|     |-∆  0|     |vn+1|
    
    
    function bem (n)
    
      x = y = zeros (n, 1) ;
    
      x(1) = 1 ; y(1) = 0 ;                       % initial conditions
      d = 2 * pi / 100 ;                          % step
    
      A = [0 d ; -d 0] ; 
      B = inv (eye (2, 2) - A) ;
    
      for t = 1 : n
        Z = B * [x(t) ; y(t)] ; x(t+1) = Z(1) ; y(t+1) = Z(2) ; 
      end
    
      plot (x, y, 'linewidth', 2)
    
    end
    

    BEM стабильнее (по крайней мере в данном случае он сходится - спираль скручивается), чем FEM, но медленнее. сумма собственных значений матрицы равна нулю (трейс), а произведение собственных значений (детерминант) равен d², значит оба собственных значения - комплексные числа (сопряженные), а значит вещественная часть равна нулю, а модуль меньше 1 - поэтому то колебания и затухают

    график на фазовой плоскости

    разница между соседними витками спирали характеризует точность

    метод трапеций

         sn+1 - sn  =  A * sn 
         
         sn+1 - sn  =  A * sn+1 
        
     теперь складываем левые и правые части и делим на два:
    
         sn+1 - sn  =  0.5 * A * (sn+1 + sn) 
    
         sn+1 - 0.5 * A * sn+1  =  sn   +   0.5 * A * sn 
    

    ЗАМЕЧАНИЕ: метод неявный, потому что искомое есть как в lhs, так и в rhs

    матрица преобразования:

    
         sn+1 = (I  -  0.5 * A)⁻ * (I  +  0.5 * A) * sn
    
    

    в октаве это может выглядеть как-то так:

    % trapezoid method
    %   u" - u = 0    
    
    %   u' =   v
    %   v' = - u
    %   |u'|     | 0    1|   |u |
    %   |  |  =  |       | * |  |
    %   |u"|     |-1    0|   |u'|  
    %
    %   un+1 = (I - (∆t/2) * A)⁻ * (I + (∆t/2) * A) * un
    
    
    function trapezoid (n)
    
      x = y = zeros (n, 1) ;
      x(1) = 1 ; y(1) = 0 ;                           % initial conditions
      d = 2 * pi / 100 ;                              % step
    
      A = [0  1 ; -1  0] ;
      I = eye (2, 2) ;
      B = (inv (I - (d / 2) * A)) * (I + (d / 2) * A) ;
    
      for t = 2 : n
        Z = B * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ;
      end
    
      plot (x, y, 'linewidth', 2) ;
    
    end
    
    матрицу нужно вычислять лишь однажды, поэтому метод быстрый. транспонированная матрица равна обратной и значит все собственные числа лежат на единичной окружности комплексной плоскости. кроме того, трейс равен нулю и значит собственные числа чисто мнимые. именно по этой причине не происходит диссипации энергии и колебания не затухают и не взлетают

    это же имеет и свои ограничения - все начинает работать не так в случае небольшой нелинейности задачи

    график на фазовой плоскости

    точность метода на каждом шаге: err ~ |u'''| * (∆t)3 / 12

    метод скачков

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

    
    пусть              u" + u  =  0  
    положим, что       u'      =  v
    тогда              v'      = -u
    
    получившаяся система:
    
            un+1 = un + ∆t * vn 
            vn+1 = vn - ∆t * un+1
    или
            un+1 = un + ∆t * vn 
            ∆t * un+1 + vn+1 = vn
    
            ( 1   0     ( un+1       ( 1  ∆t       ( un 
              ∆t  1 ) *   vn+1 )  =    0   1 )  *    vn )
    преобразуя:
    
            ( un+1       (  1     ∆t          ( un
              vn+1 )  =    -∆t  1-∆t*∆t )  *    vn )
    
    

    можно немножко покодить:

    % leap-frog method 
    %  u'' + u = 0     u(0) = 1      u'(0) = 0 
    
    %  u' = v
    %  v' = - u
    %  un+1 = un + ∆vn
    %  vn+1 = vn - ∆un+1 = vn - ∆(un + ∆vn)
    
    %  |un+1|   |  1     ∆  |   |un|
    %  |    | = |           | * |  |
    %  |vn+1|   | -∆  1-∆*∆ |   |vn|
    
    
    function leap (n)
    
      x = y = zeros (n, 1) ;
      x (1) = 1 ; y (1) = 0 ;                         % initial conditions
      d = 2 * pi / 100 ;                              % step
      A = [1 d ; -d 1-d*d] ;
    
      for t = 2 : n
        Z = A * [x(t-1) ; y(t-1)] ; x(t) = Z(1) ; y(t) = Z(2) ; 
      end
    
      plot (x, y, 'linewidth', 2) ;
    
    end
    

    график на фазовой плоскости

    Backward-Forward Difference (BFD)

      un+1 - un      un+1 - 2 * un + un-1
      ---------  +  ------------------- = A * un+1
          Δt             2 * Δt
    

    на первом шаге надо получить значение в следующей точке используя FEM. дальше - по алгоритму

    Runge-Kutta

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

    суть метода заключается в том, что раскладывают значение функции в ряд Тейлора: значение функции в нулевой точке, первая и вторая (взвешенная) производные в сумме дают достаточно точный обобщенный результат

    набросок кода выглядит приблизительно так:

    % runge-kutta second-order method
    
    %   err ~ C * ∆ ^ 3
    %   ∆ = A * x  +  A * A * x * 0.5
    
    
    %    u" = - u   
    %    u' =   v
    %    v' = - u
    
    %  |u'|     | 0  1|   |un|
    %  |  |  =  |     | * |  |
    %  |v'|     |-1  0|   |vn|
    
    
    function runge_kutta (n)
    
      x = y = zeros (n , 1) ;
      x (1) = 1 ; y (1) = 0 ;                 % initial conditions
      d = 2 * pi / 100 ;                      % step
    
      A = [0  d ; -d  0] ;
    
      for t = 2 : n
        Z0 = [x(t-1) ; y(t-1)] ;              %  the function value in x0
        Z1 = A * Z0 ;                         %  the first derivation in x0
        Z2 = A * Z1 ;                         %  the second derivation in x0
        Z  = Z0 + Z1 + 0.5 * Z2 ;             %  f(x) = f(x0) + f`(x0)*d + f"(x0)*d²/2 
        x(t) = Z(1) ; y(t) = Z(2) ;
      end
    
      plot (x, y, 'linewidth', 2) ;
    
    end
    

    результат применения этого кода:

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

    * * *

    на небольших матрицах все считается достаточно быстро. если матрицы большие, то их сначала надо преобразовать (диагонализировать, к примеру) и потом считать уже с оптимальными матрицами. но это уже другая история