Вычислить обратное преобразование Лапласа

[S_WT]

Добрый день.
По ссылке на старнице 13 есть формула (6), численное обратное преобразование Лапласа:

$ f(t) \approx \dfrac{2e^{at}}{T} \left\langle -\dfrac{1}{2} \operatorname{Re} \overline{F}(a) + \sum\limits_{k=0}^{\infty} \left\lbrace \operatorname{Re} \left\lbrace \overline{F}(a+ik\dfrac{2\pi}{T})\right\rbrace \cos k\dfrac{2\pi}{T}t - \operatorname{Im} \left\lbrace \overline{F}(a+ik\dfrac{2\pi}{T})\right\rbrace \sin k\dfrac{2\pi}{T}\right\rbrace \right\rangle $

Пытаюсь решить для момента времени t, написал следующий код:

Используется синтаксис Fortran

program inverse_laplace_transform
implicit none

!Объявление переменных
 complex, parameter     :: j = (0,1)            !Комплексная константа
 real, parameter        :: pi=3.1415926         !Число "Пи"
 complex                :: p                    !Комплексная переменная
 complex                :: w                    !Передаточная функция
 real                   :: t                    !Время
 real                   :: tau                  !Постоянная времени
 real                   :: a                    !  ???
 real                   :: T1                   !  ???
 real                   :: k, k_max              
 real                   :: lap1, lap2           
 complex                :: temp, summa, laplas

 !Начальные условия
 a = 1                  ! ???
 T1 = 1                 ! ???
 tau = 0.1              !Постоянная времени звена
 t = 0.1                !Время для которого ищем значение функции оригинала
 summa = 0                      
 k = 0
 k_max = 100000         !Большое число
 
 !Вычисление правой части формулы (суммы)      
   do while (k < k_max)
      p = a + j*k*2*pi/T1
      w = 1.0 / (tau*p + 1.0)                   
      temp = real(w) * cos(k*2*pi*t/T1) - aimag(w) * sin(k*2*pi/T1)
      summa = summa + temp
      k = k + 1
   end do
 
 !Вычисление значения преобразования Лапласа
 lap1 = 2 * exp(a*t)/T1
 p = cmplx(a, 0)
 w = 1.0 / (tau*p + 1.0)
 lap2 = -0.5 * real(w)
 laplas = lap1 * (lap2 + summa) 

print*, "laplas = ", laplas
print*, "Проверка = ", 1.0 - exp(-t/tau)

end program inverse_laplace_transform
 

Сверку делаю по передаточной функции апериодического звена

$ W(p) = \dfrac{1}{0.1p+1} $

для времени $t = 0.1$ с применением формулы

$ f(t) = 1-e^{\frac{-t}{\tau}} $

где $ \tau$ принимаю равным $0.1$

Получается следующий вывод:

 laplas   =   (1.73351192,0.00000000)
 Проверка =   0.632120550    

Результат проверки (по проверочной формуле) соответствует правде, результат обратного преобразования Лапласа нет.

Вопросы:

1. Что такое постоянная $a$, какой ее нужно принимать численно? Не понимаю что она выражает. На страницах 4-5 того же учебно-методического пособия говорится о коэффициенте $a$, но откуда он находится не указывается.
2. Что такое постоянная $T$, какой ее нужно принимать численно? То же не понимаю какой она должна приниматься и что выражает. О ней в методическом пособии говорится, что она должна быть мала, кроме того упоминается промежуток $\left\lbrace0...T\right\rbrace$.

При этом если я изменяю коэффициенты $a$ и $T$, то результаты получаются каждый раз разные. Прошу помощи.

[Farest2]

В синусе $t$ забыли (в правильности формулы не разбирался и PDF не читал).

[S_WT]

Да, спасибо, действительно просмотрел. Но исправление не помогло.
Заодно переменные temp, summa и laplas переписал в вещественный тип, хотя в данном случае это на результат не влияет.
Теперь программа выдает:

 laplas =    3.67878628    
 Проверка =   0.632120550 

[Pphantom]

Можно залезть в оригинальную статью и попробовать взять тесты из нее. Хотя в формуле из методички других опечаток, кроме уже найденной, вроде нет, но...

[S_WT]

Спасибо, к сожалению на английском, но кое что нашел. На странице 4 приводятся примеры вычисления функций и дается $T=20$. Насколько я понял это граница промежутка времени в пределах которого производится расчет точек функции $f(t)$ (их табличный расчет ведется до $T=20$). При этом указывается что $aT=5$. Пока не понятно откуда такое заключение (они пишут "мы выяснили, что для этого значения хорошо подходит $NSUM=2000$"), но принял это значение, т.е. $a = 5/20 = 0.25$. $NSUM$ принято $2000$ (у меня это k_max). В качестве функций взял их пример с подстановкой тех же самых значений:

$NSUM=2000$

$T=20$

$aT=5$ или $a=\dfrac{5}{20}=0.25$

функция:

$F(p)=p(p+1)^2$ (двойка в степени с минусом, не нашел в LaTeX помощнике как сделать степень)

$f(t) = \dfrac{t}{2}\sin(t)$

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

для $t=1$ я получаю

 laplas =   -6.45467604E-04
 Проверка =   0.420735478   

Может быть я что-то неправильно подставляю в передаточную функцию? например вот этот блок

Используется синтаксис Fortran

p = a + j*k*2*pi/T1
w = p*(p+1)**(-2.0)
temp = real(w) * cos(k*2*pi*t/T1) - aimag(w) * sin(k*2*pi*t/T1)
 

[Aritaborian]

S_WT в сообщении #1433694 писал(а):

двойка в степени с минусом, не нашел в LaTeX помощнике как сделать степень

x^{-y} $x^{-y}$.

[Farest2]

Трудно что-то посоветовать... Хоть и угадывается в исходной формуле обращение преобразования Лапласа, но в какой-то извращенной форме все записано. Так что, "прилинкованный PDF я не читал, но осуждаю". Лучше всего -- взять Диткин, Прудников "Операционное исчисление" и почитать там про обращение преобразования Лапласа: если $f(t)\Doteq F(p)$, то

$$ f(t)=\frac{1}{2\pi i}\int_{a+i{\mathbb R}} e^{tp}F(p)\,dp. $$

Интегрирование идет по вертикальной прямой, величина $a$ выбирается такой, чтобы все особые точки $F(p)$ лежали слева. То есть, для $F(p)=1/(0.1p+1)=10/(p+10)$ надо выбирать $a>-10$ (забегая вперед, лучше всего брать $a=0$, т.к. тогда $F(a)=1$). Кстати, оригинал для этого изображения вы записали неправильно (сами найдите, где ошиблись).

Переходим к интегрированию по вещественной оси:

$$ f(t)=\frac{1}{2\pi}\int_{\mathbb R} e^{t(a+ix)}F(a+ix)\,dx =\frac{e^{at}}{2\pi}\int_{\mathbb R} e^{ixt}F(a+ix)\,dx. $$

Заменяем интеграл на конечную сумму по Симпсону с шагом $h$ и получаем

$$ f(t)\approx\frac{e^{at}}{2\pi} h\sum_{-K\le k\le K} \exp(ihkt)F(a+ihk). $$

Судя по записи шага в виде $h=2\pi/T$, кто-то (авторы методички?) хотел задействовать БПФ, я буду использовать $h$ (для краткости).

Если $F(p)$ вещественна при вещественных $p$, то в сумме можно попарно собрать члены с индексами $k$ и $-k$:

$$ f(t)\approx\frac{e^{at}}{2\pi} h\Bigl\{F(a)+2\operatorname{Re}\sum_{1\le k\le K} \exp(ihkt)F(a+ihk)\Bigr\} $$

и упростить сумму (сами проделайте, никаких сопряжений у $F$ не появится).

Для вашего варианта $F(p)$ ряд при $K=\infty$ легко суммируется аналитически (ТФКП в помощь) к тому, к чему надо, а при $h\to 0$ получается и оригинал.
В общем, берите большое $K$, маленькое $h$ и все у вас получится.

[S_WT]

Farest2, спасибо за развернутй ответ.

Farest2 в сообщении #1433735 писал(а):

Интегрирование идет по вертикальной прямой, величина $a$ выбирается такой, чтобы все особые точки $F(p)$ лежали слева. То есть, для $F(p)=1/(0.1p+1)=10/(p+10)$ надо выбирать $a>-10$ (забегая вперед, лучше всего брать $a=0$, т.к. тогда $F(a)=1$).

Насколько я понимаю Вы подставляете $a$ в качестве оператора $p$ и прослеживаете, что бы функция имела область определения. В данном случае сразу возникает вопрос, в любом другой случае (имея случайно взятую передаточную функцию) $a$ нужно то же определять таким образом или можно поступать проще - брать ноль или какое-то малое положительное число?

Farest2 в сообщении #1433735 писал(а):

Кстати, оригинал для этого изображения вы записали неправильно (сами найдите, где ошиблись).

Подозреваю, что я или упростил, поскольку $k=1$, или написал знак минуса перед $t$, а не перед всей дробью, так верно будет? (здесь k это коэффициент передачи)

$f(t)=k(1-e^{-\frac{t}{\tau}})$

Farest2 в сообщении #1433735 писал(а):

Если $F(p)$ вещественна при вещественных $p$, то в сумме можно попарно собрать члены с индексами $k$ и $-k$:

Вот это, а так же на счет упрощения, мне не совсем понятно. Я к сожалению не особо математик, то что понятно математику как "самой собой" для меня может быть не понятным. Но пробую разобраться.
Пользуясь приведенной Вами формулой,

$$ f(t)\approx\frac{e^{at}}{2\pi} h\Bigl\{F(a)+2\operatorname{Re}\sum_{1\le k\le K} \exp(ihkt)F(a+ihk)\Bigr\} $$

я написал следующий код. По формуле сразу возник вопрос, $F(a)$, стоящая сразу после фигурной скобки, берется как есть или берется вещественная часть? Иначе ведь все выражение для времени $t$ получается комплексным, а так не должно быть. В коде я взял вещественную часть.

Используется синтаксис Fortran

program inverse_laplace_transform
implicit none

 complex, parameter     :: j = (0,1)                    !Комплексная константа
 real, parameter        :: pi=3.1415926                 !Число "Пи"
 complex                :: p                            !Комплексная переменная
 complex                :: w                            !Передаточная функция
 real                   :: t                            !Время
 real                   :: tau                          !Постоянная времени
 real                   :: a                            !Постоянная
 real                   :: h, k                         !Шаг и переменная k              
 real                   :: laplas
 integer                :: i, k_min, k_max
 complex                :: temp, summa

 !Исходные данные
 t = 0.1
 tau = 0.1
 a = 0
 h = 1e-3
 k_min = 1
 k_max = 100000
 summa = cmplx(0,0)
 
   !Вычисление суммы
   do i = k_min,  k_max
      k = float(i)
      p = a + j*h*k
      w = 1.0 / (tau*p + 1.0)                  
      temp = exp(j*h*k*t) * w
      summa = summa + temp
   end do
 
 !Вычисление обратного преобразования Лапласа
 p = cmplx(a, 0)
 w = 1.0 / (tau*p + 1.0)
 laplas = exp(a*t)/(2.0*pi) * h * (real (w) + 2*real(summa))

!Результат:
print'(a, f8.3)', "laplas = ", laplas
print*, "Проверка = ", 1.0 * (1.0 - exp(-t/tau))

end program inverse_laplace_transform
 

Для $t=0.1$ и при $a=0$ получаю следующее:

 laplas =    3.94315672    
 Проверка =   0.632120550  

К сожалению нет совпадения.

Кстати $a$, если его выбирать разным, все же влияет на результат. При его выборе "в пределах разумных значений" результат отличается не очень сильно, например

$a \qquad laplas $

$-5 \qquad 3.844$

$0 \qquad 3.943$

$1 \qquad 3.969 $

$5 \qquad 4.100$

Но при дальнейшем увеличении постоянный рост, например при $a=25$, $laplas = 6.229$. Поэтому все же не совсем ясно каким его нужно выбирать.
ps. Если не сочтете за труд и если это поможет, то on-line код, запускается по кнопке Run.

[Farest2]

S_WT в сообщении #1433775 писал(а):

или можно поступать проще - брать ноль или какое-то малое положительное число?

Нет, выбрать одно $a$ на все случаи жизни не получится: для $F(p)=1/(p-2)$ надо выбирать $a>2$, для $F(p)=1/(p-10^6)$ -- больше миллиона, ну а для $F(p)=1/p$ сойдет и "какое-то малое положительное число".

S_WT в сообщении #1433775 писал(а):

Подозреваю, что я или упростил, поскольку $k=1$, или написал знак минуса перед $t$, а не перед всей дробью, так верно будет?

Нет. Вы гадаете. Откройте любой учебник по операционному исчислению (того же Диткина-Прудникова) и прочитайте про оригинал изображения $F(p)=1/(p-c)$. Это минимальный порог для "вхождения" в тему, поэтому ответ я не скажу.

S_WT в сообщении #1433775 писал(а):

По формуле сразу возник вопрос, $F(a)$, стоящая сразу после фигурной скобки, берется как есть или берется вещественная часть?

Плохо дело... Ведь в вашем случае $F(p)=10/(p+10)$, т.е. для вещественных $p$ функция $F(p)$ вещественна. Зачем ей Re навешивать? Берите сразу $F(a)=10/(a+10)$, если хотите программу при разных $a$ погонять.

S_WT в сообщении #1433775 писал(а):

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

Все комплексности, которые фигурируют в правой части, съедаются Re, которая стоит перед суммой. В итоге остается только вещественное выражение. Про слагаемое $F(a)$ я уже сказал. Как упростить сумму -- это сами. Ну, ведь знаете же вы, что $\exp(ix)=\cos(x)+i\sin(x)$, или нет? И Re- и Im-части функции $F(a+ihk)=10/(a+ihk+10)$ сами должны найти.

Программу даже смотреть не буду. Пока вы не поймете математику, которую собираетесь перевести в код, нет смысла. А когда поймете, то и я уже буду не нужен.

[S_WT]

Я спутал весовую и переходную функции. Все работает.

При $kmax=1000000$ и для $\tau = 0.1$ для переходной функции получились следующие результаты:

t            laplas      точное знач.   разница по модулю
==========================================================
0.0010      0.0097      0.0100         2.765E-04
0.0100      0.0944      0.0952         7.849E-04
0.0500      0.3932      0.3935         2.234E-04
0.1000      0.6321      0.6321         1.591E-05
0.2000      0.8644      0.8647         2.689E-04
0.3000      0.9503      0.9502         7.004E-05
0.4000      0.9816      0.9817         5.347E-05
0.5000      0.9933      0.9933         4.572E-05
1.0000      0.9999      1.0000         3.988E-05

Можно попробовать поэкспериментировать с уменьшением шага, при расчете приведенных данных он был принят $1e-3$. Зависимость от коэффициента $a$ сильно уменьшилась (возможно потому, что в данном случае производится расчет переходной функции и разброс значений стал менее заметен), однако при больших значениях все равно погрешность возрастает, насколько я понял его нужно выбирать весьма умеренно. Можно и $kmax$ попробовать задрать еще выше.

ps Я это проходил около 30 лет назад и с тех пор понадобилось в первый раз, и то не по работе.

[Farest2]

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

S_WT в сообщении #1433815 писал(а):

при больших значениях все равно погрешность возрастает

то это вполне предсказуемо: $f(t)$ -- хорошая убывающая функция, а ей в формуле множитель $e^{at}$ подсовывают. Ясно, что при больших $a$ и $t$ функция $f(t)$ вынуждена вычисляться как что-то типа $\infty\cdot 0$, что не способствует точности вычислений.

[novichok2018]

Работающие обращения пр. Лапласа - это достаточно большая и разработанная область, намного шире, чем указано в методичке, всего пара методов. Если интересно с этим разобраться, то можно начать с обращения чисто действительного в книге Рамма "Многомерные обратные задачи рассеяния". Из другого посмотреть инет и например работы Saito по теме.

[S_WT]

Farest2 в сообщении #1433979 писал(а):

то это вполне предсказуемо: $f(t)$ -- хорошая убывающая функция, а ей в формуле множитель $e^{at}$ подсовывают. Ясно, что при больших $a$ и $t$ функция $f(t)$ вынуждена вычисляться как что-то типа $\infty\cdot 0$, что не способствует точности вычислений.

Можно включить двойную точность еще. Тем не менее протестировал на примерах из интернета с разными передаточными функциями, графики переходных процессов везде получились полностью идентичными, $kmax$ было установлено 500 000. Любопытно, каким устанавливается этот параметр в других реализациях…

novichok2018 в сообщении #1434036 писал(а):

Работающие обращения пр. Лапласа - это достаточно большая и разработанная область, намного шире, чем указано в методичке,

Я видел в методичке, хотя как Вы говорите там вероятно не все, выбрал наиболее простой способ. Читал, что можно и с помощью преобразования Фурье отыскивать значение оригинала. Какие у него плюсы и минусы?

[S_WT]

Провел анализ требуемого максимального $k$ посредством построеня графиков переходных процессов для колебательного звена. Были построены графики при $k=10000$, $50000$, $100000$, $250000$, $500000$ и $1000000$. При $k=10000$ искажения графика очень сильные (вплоть до неузнаваемости), при $k=50000$ искажений много меньше, но они достаточно сильные. При $k=100000$ график принимает практически идеальную форму, но все же есть его незначительные отклонения, при $k=250000$, $500000$ и $1000000$ графики не отличимы друг от друга. Иначе говоря сделан вывод, что достаточно применять $k=250000$. Только для построения некоторых графиков (например имеющих какие-то сильные колебания, выбросы и т.п.) вероятно имеет смысл, на всякий случай, принимать $k$ больший, чем 250000. При построении всех графиков было принято значение $h=0.001$. Вероятно есть смысл попробовать и с меньшим значением $h$.