ラプラス変換で<1/3>微分方程式 「慶應大学講義 物理情報数学C 第十回 ラプラス変換の応用,信号...」をChatGPTとWolframAlphaとsympyでやってみたい。

・微分方程式をラプラス変換で。

・別頁<２/3>,<3/3>と続きます。
・2023/10/18作成中

オリジナル

(0:00〜6:06) 微分方程式

< Modeling and control of dynamic systems - From electric vehicles to systems biology : Adachi's Group
https://youtu.be/ZEWz1mMVJgY

(00:00) ＞時間領域。
(00:00) ＞ラプラス領域。
(04:51) ＞これ答え。

ChatGPT-3.5先生へ。(???わかりませんでした。???)

・???微分法方程式。中国語になりました。小技かも。???

微分方程式
x''(t)+5*x'(t)+4*x(t)=0
x(0)=2,x'(0)=-1

かっこ(t)なし。

微分方程式
x''+5*x'+4*x=0
x(0)=2,x'(0)=-1

(省略)

Mathematicaで

いつの日か、実行してみたい。

WolframAlphaで(できました。)

上と同じ。かっこ(t)なし。

sympyで

・???　tは,nonnegativeですか　？？？
・sympyの自作myRld関数(小数を分数変換)で、２回実行した。

# ver0.1
from sympy import *
import re
var('s',real    =True)
var('t',positive=True)
x,X=map(Function,('x','X'))
def myRld(f): # ver0.1
     li = re.findall(r'[.0-9]+',str(f))
     sorted_li=sorted(li ,key=lambda s:len(s))
     syosu=                 sorted_li[-1]
     bunsu=Rational(sympify(syosu)).limit_denominator(10**12)
     return sympify( str(f).replace(str(syosu),str(bunsu)) )
def myLapDiff(xdash,x0,x10):
     if   xdash==2: return s**2*X(s)-s*x0-x10
     elif xdash==1: return s**1*X(s)-  x0
     elif xdash==0: return      X(s)
     else:          return print("#error myLapDiff")
def myLap(g):
     return laplace_transform(g,t,s)[0]
def myInvLap(g):
     return inverse_laplace_transform(g,s,t).evalf().simplify()
#---------------------------------------------------------------------------------------------------------------------
eq01=Eq(     diff(x(t),t,2    )+5     *diff(x(t),t,1    )+4     *diff(x(t),t,0  ),     0  ) ;print("#(00:00)01",eq01)
print()
eq11=Eq(myLapDiff(       2,2,-1)+5*myLapDiff(     1,2,-1)+4*myLapDiff(    0,2,-1),myLap(0)) ;print("#(00:00)11",eq11)
eq12=Eq(X(s),      solve(eq11,X(s))[0])                                                     ;print("#(03:00)12",eq12)
eq13=Eq(X(s),      apart(eq12.rhs)    )                                                     ;print("#(03:11)13",eq13)
eq14=Eq(x(t),   myInvLap(eq13.rhs)    )                                                     ;print("#(04:36)14",eq14)
eq15=Eq(x(t),myRld(myRld(eq14.rhs)    ))                                                    ;print("#(04:51)15",eq15)
#(00:00)01 Eq(4*x(t) + 5*Derivative(x(t), t) + Derivative(x(t), (t, 2)), 0)

#(00:00)11 Eq(s**2*X(s) + 5*s*X(s) - 2*s + 4*X(s) - 9, 0)
#(03:00)12 Eq(X(s), (2*s + 9)/(s**2 + 5*s + 4))
#(03:11)13 Eq(X(s), -1/(3*(s + 4)) + 7/(3*(s + 1)))
#(04:36)14 Eq(x(t), 2.33333333333333*exp(-t) - 0.333333333333333*exp(-4*t))
#(04:51)15 Eq(x(t), 7*exp(-t)/3 - exp(-4*t)/3)

いつもの? sympyの実行環境と 参考のおすすめです。

(テンプレート)

参考

参考文献

材料力学演習
小井土正六著. -- 学献社, 1963.
5.4 ラプラス変換を利用する法
P138
https://opactmu.lib.tmu.ac.jp/iwjs0013vir/vhol/holdetail/?selectedHolid=HL01134776

次頁<２/3>単振動へ続く