ラプラス変換で<2/3>単振動:力学システムの例 「慶應大学講義 物理情報数学C 第十回 ラプラス変換の応用,信号...」をChatGPTとWolframAlphaとsympyでやってみたい。

前頁<1/3>微分方程式より

・単振動をラプラス変換で。

オリジナル

(6:06〜13:40) 単振動

(07:50) ＞静かに離す。速度の初期値が0。
(10:28) ＞因数分解したくない.オイラーの公式でできる。 ???
(10:00) ＞絶対暗記ね。
(11:12) ＞固有周波数(英:natural frequency)
(13:07) ＞単振動

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

微分方程式
m*x''(t)+k*x(t)=0,x(0)=0,x'(0)=0

かっこ(t)なし。

微分方程式
m*x''+k*x=0 ,x(0)=0,x'(0)=0

(省略)

Mathematicaで

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

WolframAlphaで(???わかりませんでした。???)

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

標準の計算時間制限を超えました...

sympyで

・ver0.13 三角関数は、暗記しなくてもいいですか？

# ver0.13
from sympy import *
var('s'  ,real    =True)
var('t'  ,positive=True)
var('m,k',real    =True)
var('x0' ,real    =True)
var('w0' ,real    =True)
x,X=map(Function,('x','X'))
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()
def myY2y(Ys):
     return myInvLap(apart(Ys))
# 絶対暗記ね。 (10:00)
repw={sqrt(k/m)**2  :w0**2    }  
# repc={s/(s**2+w0**2):cos(w0*t)}
#------------------------------------------------------------------------------------------------- 
eq01=Eq(m*     diff(x(t),t,2     )+k     *diff(x(t),t,0     ),      0  ) ;print("#(08:00)01",eq01)
print()
eq11=Eq(m*myLapDiff(       2,x0,0)+k*myLapDiff(       0,x0,0),myLap(0) ) ;print("#(08:00)11",eq11)
eq12=Eq(X(s),solve(eq11,X(s))[0] )                                       ;print("#(08:46)12",eq12)
eq13=Eq(eq12.lhs,(numer(eq12.rhs)/m) /( (denom(eq12.rhs)/m).expand() ) ) ;print("#(09:24)13",eq13)
eq14=Eq(eq13.lhs,eq13.rhs.subs(repw))                                    ;print("#(09:24)14",eq14)
eq15=Eq(x(t)    ,myInvLap(eq14.rhs ))                                    ;print("#(09:24)15",eq15)
# 途中省略
#(09:24)15 Eq(x(t), x0*cos(t*w0))

・ver0.12

# ver0.12
from sympy import *
var('s'  ,real    =True)
var('t'  ,positive=True)
var('m,k',real    =True)
var('x0' ,real    =True)
var('w0' ,real    =True)
x,X=map(Function,('x','X'))
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()
def myY2y(Ys):
     return myInvLap(apart(Ys))
# 絶対暗記ね。 (10:00)
repw={sqrt(k/m)**2  :w0**2    }  
repc={s/(s**2+w0**2):cos(w0*t)}
#------------------------------------------------------------------------------------------------- 
eq01=Eq(m*     diff(x(t),t,2     )+k     *diff(x(t),t,0     ),      0  ) ;print("#(08:00)01",eq01)
print()
eq11=Eq(m*myLapDiff(       2,x0,0)+k*myLapDiff(       0,x0,0),myLap(0) ) ;print("#(08:00)11",eq11)
eq12=Eq(X(s),solve(eq11,X(s))[0] )                                       ;print("#(08:46)12",eq12)
eq13=Eq(eq12.lhs,(numer(eq12.rhs)/m) /( (denom(eq12.rhs)/m).expand() ) ) ;print("#(09:24)13",eq13)
eq14=Eq(eq13.lhs,eq13.rhs.subs(repw))                                    ;print("#(09:24)14",eq14)
eq15=Eq(x(t)    ,eq14.rhs.subs(repc))                                    ;print("#(09:24)15",eq15)
#(08:00)01 Eq(k*x(t) + m*Derivative(x(t), (t, 2)), 0)

#(08:00)11 Eq(k*X(s) + m*(s**2*X(s) - s*x0), 0)
#(08:46)12 Eq(X(s), m*s*x0/(k + m*s**2))
#(09:24)13 Eq(X(s), s*x0/(k/m + s**2))
#(09:24)14 Eq(X(s), s*x0/(s**2 + w0**2))
#(09:24)15 Eq(x(t), x0*cos(t*w0))

・ver0.11
・終了しません。全部暗記しないは、ダメでした。

# ver0.11
from sympy import *
var('s'  ,real    =True)
var('t'  ,positive=True)
var('m,k',real    =True)
var('x0' ,real    =True)
var('w0' ,real    =True)
x,X=map(Function,('x','X'))
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()
def myY2y(Ys):
     return myInvLap(apart(Ys))
# 絶対暗記ね。 (10:00)
# repw={sqrt(k/m)**2  :w0**2    }  
# repc={s/(s**2+w0**2):cos(w0*t)}
#------------------------------------------------------------------------------------------------- 
eq01=Eq(m*     diff(x(t),t,2     )+k     *diff(x(t),t,0     ),      0  ) ;print("#(08:00)01",eq01)
print()
eq11=Eq(m*myLapDiff(       2,x0,0)+k*myLapDiff(       0,x0,0),myLap(0) ) ;print("#(08:00)11",eq11)
eq12=Eq(X(s),solve(eq11,X(s))[0] )                                       ;print("#(08:46)12",eq12)
# eq13=Eq(eq12.lhs,(numer(eq12.rhs)/m) /( (denom(eq12.rhs)/m).expand() ) ) ;print("#(09:24)13",eq13)
# eq14=Eq(eq13.lhs,eq13.rhs.subs(repw))                                    ;print("#(09:24)14",eq14)
# eq15A=Eq(x(t)    ,      myInvLap(eq12.rhs) )                             ;print("#(09:24)15",eq15A)
eq15B=Eq(x(t)    ,myY2y(myInvLap(eq12.rhs)))                             ;print("#(09:24)15",eq15B)
# eq15A、eq15A
# 終了しません。

[次頁<3/3>電気回路へ続く]作業中