сигналы     преобразование Фурье

преобразование Лапласа


  • правило экспоненциального сдвига
  • правило подобия
  • правило получения образа производной

  • как все начиналось?

               ∞
               Σ aₙ * xⁿ   =  C (x)
              n=0
    
               ∞
               Σ a (n) * xⁿ  = C (x)
              n=0
    
            ----------------------------------------------
             дискр ф-ция                  непрер ф-ция
            ----------------------------------------------
             a(n)           ~~~~~~~》          C(x)
            ----------------------------------------------
              1             ~~~~~~~》   1/(1-x) , |x| < 1
    
             1/n!           ~~~~~~~》        exp (x)
            ----------------------------------------------      
    
             ∞
             ∫ f(t) * x^t dt   =  F (x) , 0 < x < 1 
             0
    
             ∞
             ∫ f(t) * [exp (ln x)]^t dt   =  F (x)  ,  ln x < 0
             0
    
             ∞
             ∫ f(t) * exp (-s * t) dt   =  F (s)  ,  s = - ln x 
             0 


    линейный оператор, действующий в пространстве функций f от вещественного положительного аргумента t:

                             ∞
                   ƚ(s)  =   ∫ f(t) * exp(-s * t)  dt
                             0
    называется "преобразованием Лапласа". он отображает множество функций неотрицательной вещественной переменной t на множество функций (вообще говоря - комплексной) переменной s

    ограничения - функция f(t) должна быть дифференцируема и интеграл должен сходиться. "плохими" функциями являются, например, функции 1/t и exp(t²)

    посчитаем для простых функций, для которых есть аналитическое решение:

        тривиальный случай: f = 0
    
            ∞
        ƚ = ∫ 0 * exp(-s*t) dt = 0
            0
    
    
    обобщенная функция Дирака: δ₀ очевидно, что ƚ δ₀ = 1 ибо ∞ +0 +0 ∫ δ₀ * exp (-s*t) dt = ∫ δ₀ * exp (-s*0) dt = ∫ δ₀ dt = 1 0 -0 -0
    функция Хевисайда: f = 1 ∞ ∞ ƚ = ∫ 1 * exp (-s*t) dt = - exp(-s*t) / s | = 1/s 0 0
    функция склона: f = t ∞ ∞ ∞ ƚ = ∫ t * exp (-s*t) dt = [- exp(-s*t) / s ] * t | - ∫ 1 * [- exp (-s*t) / s] dt 0 0 0 ∞ ƚ = 0 - [1/s * - exp (-s*t)/s] | = 1/s² 0
    наша любимая экспонента: f = exp(t) ∞ ∞ ƚ = ∫ exp (-s*t) * exp(t) dt = - exp(-(s-1)*t) / (s - 1) | = 1 / (s - 1) 0 0


    в СКА Максима:


    (%i86) laplace (exp(t),t,s);
                                           1
    (%o86)                               -----
                                         s - 1
    (%i87) laplace (t,t,s);
                                          1
    (%o87)                                --
                                           2
                                          s
    (%i88) laplace (sin(t),t,s);
                                          1
    (%o88)                              ------
                                         2
                                        s  + 1
    (%i89) laplace (cos(t),t,s);
                                          s
    (%o89)                              ------
                                         2
                                        s  + 1
    (%i90) laplace (exp(-t),t,s);
                                           1
    (%o90)                               -----
                                         s + 1
    (%i91) laplace (-t,t,s);
                                           1
    (%o91)                               - --
                                            2
                                           s
    
    (%i92) ilt (1/s^2,s,t);
    (%o92)                                 t
    (%i93) ilt (1/(s-1),s,t);
                                            t
    (%o93)                                %e
    


    правило экспоненциального сдвига

                                    ∞
           ℒ [exp (a * t) * f(t)] = ∫ f(t) * exp [- (s - a) * t] dt = ƚ (s - a)
                                    0
    
                где ƚ(s) есть ℒ [f (t)]

    но тогда сразу же, полагая f (t) = 1, получим:

          ℒ [cos(t)] = ℒ [(exp(i*t) * 1 + exp(-i*t) * 1) / 2]
    
                     = [1 / (s - i) + 1 / (s + i)] / 2
    
                     = s / (s² + 1)
    
    и аналогично
    
          ℒ [sin(t)] = ℒ [(exp(i*t) * 1 - exp(-i*t) * 1) / 2i]
    
                     = [1 / (s - i) - 1 / (s + i)] / 2i
    
                     = 1 / (s² + 1)


    правило подобия

           ℒ [f (k * t)] = 1/k * ƚ (s / k)
    
           где ƚ (s) есть ℒ [f(t)],  k ∈ ℝ₊

    действительно

                         ∞
           ℒ [f(k*t))] = ∫ f(k*t) * exp(-s*t) dt =
                         0
    
           пусть k*t = u
    
                                ∞
                       =  1/k * ∫ exp(-s*u/k) * f(u) du
                                0

    но тогда сразу же

             exp (k * t) = 1/k * [ 1  / (s/k - 1) ]        =  1 / (s  - k)
    
             cos (k * t) = 1/k * [ (s/k) / ((s/k)² + 1) ]  =  s / (s² + k²)
    
             sin (k * t) = 1/k * [    1  / ((s/k)² + 1) ]  =  k / (s² + k²)

    таблица образов

    обратное преобразование Лапласа вычисляется по формуле:

               ∞
        f(t) = ∫ F(s) * exp (s * t) ds
              -∞  

    но всегда для получения прообраза из образа пользуются таблицами:
    оригинал : f
    временной домен
    образ : s
    частотный домен
    обобщенные
    0 0
    δ₀ 1
    δ₁ 1/s
    t 1/s²
    экспоненты
    exp (k * t) 1 / (s - k)
    t * exp (k * t) 1 / (s - k)²
    тригонометрические
    cos (k * t) s / (s² + k²)
    sin (k * t) k / (s² + k²)
    exp (a * t) * cos (k * t) (s - a) / [ (s - a)² + k²]
    exp (a * t) * sin (k * t) k / [ (s - a)² + k²]

    этого вполне достаточно для инженерных расчетов цепей


    теорема об образе производной

            ℒ [ f'(t) ]   =   s * ƚ(f)  -  f(0)
    
            ℒ [ f˝(t) ]   =   s² * ƚ(f)  -  s * f'(0)  -  f(0)


    примеры

    пример расчета RC-цепи

    есть сигнал на входе, требуется вычислить сигнал на выходе

    
         uc    +    ur            = v(t)
         uc    +    ic * R        = v(t)
         uc    +    u'c * C * R   = v(t)
    
    
                u'        +    u / R*C     =    v(t) / R*C
    пусть
         u(0)   =  0
         R * C  =  1
         v(t)   =  exp (-2 * t)
    тогда
         s *  ƚ(f)   -   0    +   ƚ(f)   =   1 / (s + 2)
    
         ƚ(f) * (s + 1)                  =   1 / (s + 2)
    
         ƚ(f)                            =   1 / [(s + 1) * (s + 2)]
    
    
         первый полюс: -1
         второй полюс: -2
    
         поскольку степень числителя меньше степени знаменателя:
         производная знаменателя (s² + 3 * s + 2)' = 2 * s + 3
         в первом полюсе равна:   1
         во втором полюсе равна: -1
    
         результат:
    
                   1        -1
         ƚ(f) =  -----  +  -----
                 s + 1     s + 2
    
    ответ: u(t) = exp(-1 * t) - exp(-2 * t)

    listing octave:


    function y = rccirc (t)
      y = exp(-1 * t) - exp(-2 * t);
    end
    
    t = linspace (0, 6, 1000)
    
    plot (t, rccirc(t), 'linewidth', 2, t, exp (-2*t), 'linewidth', 2)
      

    результат:


    пример расчета RL-цепи

    есть сигнал на входе, требуется вычислить сигнал на выходе

         ul       +   ur          = v(t)
         L * i'l  +   il * R      = v(t)
         i'l      +   il * R/L    = v(t) / L
    
    
            i'        +      i * R/L       =      v(t) / L
    
    пусть
         i(0)   =  0
         R      =  2
         L      =  1
         v(t)   =  3
    тогда
         s * ƚ(f)  -   0   +   2 * ƚ(f)  =   3 / s
    
         ƚ(f)  *  (s + 2)                =   3 / s
    
         ƚ(f)                            =   3 / [s * (s + 2)]
    
         первый полюс:  0
         второй полюс: -2
    
         поскольку степень числителя меньше степени знаменателя:
         производная знаменателя (s² + 2 * s)' = 2 * s + 2
         в первом полюсе равна:    2
         во втором полюсе равна:  -2
    
         результат:
    
                 1.5     1.5
         ƚ(f) =  ---  - -----
                  s     s + 2
    
    ответ: i(t) = 1.5 * (1 - exp(-2 * t))

    listing octave:


    function y = rlcirc (t)
      y = 1.5 * (1 - exp (-2 * t));
    end
    
    t = linspace (0, 6, 1000)
    
    plot (t, rlcirc(t), 'linewidth', 2, [0,6], [3,3], 'linewidth', 2)
      

    результат:


    соотношения для элементов цепей
     для резистора:      V(s) = R * I(s)
    
     для конденсатора:   V(s) = I(s) / (C*s)  +  v₀ / s
                         I(s) = (C*s) * V(s)  -  C * v₀
                         емкость предсталена элементом с импедансом 1/Cs
                         последовательно к которому подключен источник напряжения v₀/s
    
     для катушки:        I(s) = V(s) / (L*s)  +  i₀ / s
                         V(s) = (L*s) * I(s)  -  L * i₀
                         индуктивность представлена элементом с импедансом Ls
                         последовательно к которому подключен источник тока i₀/s
    
     где                 V(s) = ƚ[u(t)]
                         I(s) = ƚ[i(t)]
                         v₀, i₀ - начальные условия
    

    дифуры второго порядка

    теперь, поскольку преобразование Лапласа - линейный оператор, то для постоянно встречающихся в инженерных расчетах уравнений второго порядка имеем соответствие:

                f˝(t)                 +        k₁ * f'(t)          +    k₂ * f(t)
    
                 ↓                              ↓                       ↓
    
       s² * ƚ(f) - s * f(0) - f'(0)   +    k₁ * s * ƚ(f) - f(0)    +    k₂ * ƚ(f)
    
    
         f˝(t) + k₁ * f'(t) + k₂ * f(t)      →        (s² + k₁*s + k₂) * ƚ (f)
                                                  -            (s - 1) * f (0)
                                                  -                      f'(0)
    

    аналитический пример


    пример расчета параллельного контура

    пусть R = 1 L = 1 C = 1 v₀ = 0 i₀ = 0 v = δ₀ , тогда

                   V(s) = 1
    
                                  1
                   I(s) = -------------------
                          R + (L*s) + 1/(C*s)
    
                                   Cs
                   I(s) = --------------------
                          (C*R)s + (L*C)s² + 1
    тогда
    
                               s                          s
                   I(s) = ----------- = --------------------------------------
                          s² + s + 1    (s - (-.5 + .87i) * (s - (-.5 - .87i))
    
                           -1.5 - .87i         -1.5 + .87i
                   I(s) = --------------  +  ---------------
                          s + (.5 -.87i)     s + (.5 + .87i)
    
        Ответ:
    
                   i(t) = (-1.5 - .87i) * exp[(-.5 + .87i)*t] + (-1.5 + .87i) * exp[(-.5 - .87i)*t]
    

    закодим это дело, listing octave:


    function y = rlcp (t)
      y = (-1.5 - .87i) * exp((-.5 + .87i)*t) + (-1.5 + .87i) * exp((-.5 - .87i)*t);
    end
    
    t = linspace (0, 30, 1000);
    
    plot (t, rlcp(t), "r", 'linewidth', 2)
      

    получим:

    лепота!

    положим теперь, что v = δ₁

        тогда      V(s) = s
    
                               s            s²                         s²
                   I(s) = ----------- = ---------- = --------------------------------------
                          1 + s + 1/s   s² + s + 1   (s - (-.5 + .87i) * (s - (-.5 - .87i))
    
                             -.5 - .3i            -.5 + .3i
                   I(s) = ----------------- + ------------------
                          (s - (-.5 + .87i)   (s - (-.5 - .87i))
    
        ответ:
    
                   i(t) = (-.5 - .3i) * exp[(-.5 + .87i) * t] + (-0.5 + .3i) * exp[(-.5 -.87i) * t]
    

    закодим это дело. listing octave:


    function y = rlcp (t)
      y = (-.5 - .3i) * exp((-.5 + .87i)*t) + (-0.5 + .3i) * exp((-.5 - .87i)*t);
    end
    
    t = linspace (0, 30, 1000);
    
    plot (t, rlcp(t), "r", 'linewidth', 2)
      

    получим:

    красота!

    положим теперь, что v = δ₀ , а R = 0

    тогда  V(s) = 1
    
                            1          s
                   I(s) = ------- = ------
                          s + 1/s   s² + 1
    
        ответ:
    
                   i(t) = cos(t)
    

    положим теперь, что v = δ₁ , при R = 0

    тогда  V(s) = s
    
                             s         s²
                   I(s) = ------- = ------ ≈ 1
                          s + 1/s   s² + 1
    
        ответ:
    
                   i(t) ≈ δ₀
    


    пример расчета последовательного контура

    TODO
    

    пример расчета последовательно-параллельного контура

    TODO
    

    итак, мы используем преобразование Лапласа для решения неоднородных дифференциальных уравнений. записываем в левой части производные и функции, а в правой части - свободные члены:

        y˝(t) + y'(t) + y(t) = f(t)
        y˝(0) = a2
        y'(0) = a1
        y (0) = a0
      
    затем делаем преобразования Лапласа для левой и правой части (учитывая начальные условия), находим представление дробно-рационального образа в виде суммы простейших и производим обратное преобразование Лапласа (обычно из таблицы образов) - получая тем самым y(t)

    преобразование Лапласса и автономные линейные системы

    что есть преобразование Лапласса? это линейный оператор над пространством всюду дифференцируемых функций от неотрицательного вещественного числа

    пусть есть функция f. линейный оператор ℒ - это оператор Лапласса и

           s  =  ℒ (f)
    s называется образом функции f

    всегда выполняется равенство:

        ℒ (f') = s * ℒ (f)  -  f (0)  

    обратное преобразование Лапласса переводит образ функции обратно в функцию:

        ℒ⁻ (s) = ℒ⁻ (ℒ (f)) = f 

    в канонической форме:
    x - вектор состояний системы
    A - матрица автономной системы:

        x'  =  A * x

    применяем преобразование Лапласса к левой и правой части:

                 s * xs - x₀  =  A * xs
    
              (s * I - A) * xs  =  x₀
    
                            xs  =  (s * I  -  A)⁻ * x₀

    матрица (s * I - A)⁻ может быть получена, если det (A) ≠ 0

    теперь применим к левой и правой частям обратное преобразование Лапласса. оно тоже линейно и поэтому коэфициент x₀ просто выносится за знак операции:

                            xt  =   ℒ⁻ (s * I - A)⁻  *  x₀
    пусть
              Фt  =   ℒ⁻ (s * I - A)⁻  
    тогда
              xt  =   Фt * x₀

    матрица Ф(t) = ℒ⁻ X(s) называется матрица перехода

    преобразования образов в оригиналы берутся из таблиц и производятся для матрицы (s * I - A)⁻ поэлементно

    пример c интегратором

    решим численно такую систему:

         x₁' = x₂
         x₂' = 0

    в канонической форме:

      ( x₁'       ( 0 1      ( x₁
        x₂' )  =    0 0 ) *    x₂ )
    
      s * I - A  = ( s  -1
                     0   s )
    
      (s * I - A)⁻  = ( 1/s   1/s²
                         0    1/s  )

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

      Ф(t)  =  ℒ⁻ ( 1/s   1/s²
                     0    1/s  )
    
      Ф(t)  =  ( 1   t
                 0   1 )
    следовательно, решение:
        x(t)  =  [1  t; 0  1] * x0 

    кодируем, выбираем начальные условия и запускаемся. все просто

    пример с осцилятором

    решим численно такую вот систему:

         x₁' =  x₂
         x₂' = -x₁ 

    в канонической форме:

      ( x₁'      ( 0 1     ( x₁
        x₂' ) =   -1 0 ) *   x₂ )
    
      s * I - A = ( s -1
                    1  s )
    
      (s* I  -  A)⁻ = ( s/(s²+1)   1/(s²+1)
                       -1/(s²+1)   s/(s²+1) )
    
    
      Ф(t)  =  ℒ⁻ ( s/(s²+1)  1/(s²+1)
                   -1/(s²+1)  s/(s²+1) )
    
      Ф(t)  =  ( cos(t)  sin(t)
                -sin(t)  cos(t) )
    следовательно, решение:
        x(t)  =  [cos(t)  sin(t) ; -sin(t) cos(t)] * x0

    кодируем, выбираем начальные условия, запускаемся, получаем график x(t)

    а теперь - самое интересное...

    экспонента матрицы

              (I - A)⁻ = I  +  A  +  A * A  +  A * A * A  + ..
    

    если собственные значения А меньше 1, то такой ряд сходится

    теперь

              (s * I  -  A)⁻  =  1/s * (I  -  A ./s)⁻
    
                              =  I ./s  +  A ./s²  +  A*A ./s³  + ...
    
    
              Ф(t) = ℒ⁻ (I ./s  +   A ./s²  +   A*A ./s³      + ... )
    
              Ф(t) =     I      +   t * A   +  (t * A)² ./2!  +  (t * A)³ ./3! + ...
    
              Ф(t) = expm (t * A) 

    это сильно упрощает дело кодирования. ок, ну решим-ка теперь предыдущие два примера численно. сначала накодим:

    %% inputs:
    %%    A  - nxn matrix of LDS
    %%    x0 - nx1 vector of initial state
    %%    t  - nx1 vector of time steps
    %% output
    %%       - plot x'(t)     each component
    
    function exponm100 (A, x0, t)
      for i = 1 : columns (t)
        y(i, :) = expm (A * t(i)) * x0;
      end
      plot (t, y, 'linewidth', 2);
    end

    запускаем с матрицей интегратора

    octave> t = linspace (0,10,100);
    octave> A = [0 1; 0 0]
    A =
    
       0   1
       0   0
    
    octave> x0
    x0 =
    
       1
       1
    
    octave> exponm100 (A, x0, t)
    получаем результат

    запускаем с матрицей осцилятора

    octave> A = [0  1; -1  0]
    A =
    
       0   1
      -1   0
    
    octave> x0 = [ 1; 1]
    x0 =
    
       1
       1
    
    octave> t = linspace (0,10,100);
    octave> exponm100 (A, x0, t)

    получаем результат

    забываем о таблицах образов функций, как о страшном сне, и помним, что ныне и приисно и во веки веков:

                              x'  =  A * x
    
                              xt  =  expm (t * A) * x₀
    

    теперь, очень важное замечание!

                expm (A + B)  ≠ (expm A) * (expm B)
    
                но, если A * B = B * A
                то expm (A + B) = (expm A) * (expm B)
    
    
                и, значит, expm (t₁ * A  +  t₂ * A) = expm (t₁ * A)  *  expm (t₂ * A)

    применим это:

            x(τ + t) = expm ((t + τ) * A)
    
            x(τ + t) = expm (t * A  + τ * A)
    
            x(τ + t) = expm (t * A) * expm (τ * A)
    
            x(τ + t) = expm (t * A) * x(τ)

    значит, если мы знаем состояние системы в какой-то момент времени τ мы сможем вычислить состояние системы через t тиков или же мы можем вычислить, каким было состояние системы t тиков назад (взяв t с минусом)

    давайте проверимся на предыдущем примере интегратора

    в момент времени τ = 8 состояние системы было x(8) = [9; 1]. значит для t = -8 должно получиться х(0) = [1; 1]

    проверим, так ли это:

    octave> A = [0 1; 0 0];
    octave> x8 = [9; 1];
    octave> t8 = linspace (0, -8, 80);
    octave> exponm100 (A, x8, t8)
    

    получаем результат

    давайте теперь проверимся на предыдущем примере осцилатора

    для осцилятора x(8) = [0.85; -1.40]. для t = -8 должно получиться х(0) = [1; 1]

    проверим

    octave> A = [0 1; -1 0];
    octave> x8 = [ -1.145; 0.843];
    octave> t8 = linspace (8,0,80);
    octave> exponm100 (A, x8, t8)
    

    получаем результат

    шайтан-машина времени, однако получилась!

    еще один пример

    итак, если есть система дифур:

                v'    =    A * v
    при условии, что А - квадратная, полного ранга:
              v(t)    =    expm (t * A)   *    v₀

    дана система:

                ( x'      (  -2.00    0.25    0.10       ( x
    
                  y'          9.00   -2.00    0.00         y
    
                  z' ) =     -1.00    3.00   -4.00 )  *    z )

    итак:

    function firstODE (A, v0, t, m)
      n = rows (A);
      if (n != rank(A)) disp("det A iz zero - quiting"); exit; end
    
      t = 0 : t/m : t;
      y = zeros (n, m);
      y(:, 1) = v0;
    
      for i = 2 : m
        y(:, i) = expm (t(i) * A) * v0;
      end
    
      for i = 1 : n
        subplot (3, 1, i)
        plot (x = 0 : m - 1, y (i, x + 1), 'linewidth', 2);
        grid on
      end
    end

    запускаем код:

    octave> A
    A =
    
      -2.00000   0.25000   0.10000
       9.00000  -2.00000   0.00000
      -1.00000   3.00000  -4.00000
    
    octave> v0
    v0 =
    
       1
       0
       0
    
    octave> firstODE (A, v0, 7, 100)

    получаем результат: