為了更清楚，我想繪制鐘擺阻尼振蕩的二階微分方程的解。鏈接到有關(guān)所用方程式的維基百科：https ://en.wikipedia.org/wiki/Harmonic_oscillatorfrom sympy.interactive import printingprinting.init_printing(use_latex=True)import numpy as npimport scipy as spfrom sympy import*mport sympy as syp`from scipy.integrate import odeintimport matplotlib.pyplot as pltt=syp.Symbol('t')x=syp.Function('x')(t)m=2.0k=5.0a=0.5z=a/(2.0*np.sqrt(m*k))w=np.sqrt(k/m)eq=x.diff(t,t)+2.0*z*w*x.diff(t)+w**2.0*xdsolve(eq,t,0,ics={eq(1.0):0,eq(2.0):5})
1 回答
慕哥9229398
TA貢獻(xiàn)1877條經(jīng)驗(yàn) 獲得超6個(gè)贊
您沒有ics按預(yù)期構(gòu)造參數(shù)：
In [6]: dsolve(eq, ics={x.subs(t, 1.0): 0, x.subs(t, 2.0): 5})
Out[6]:
-0.125?t
x(t) = (-0.0346285740992263?sin(1.57619002661481?t) - 6.42012708343871?cos(1.57619002661481?t))??
如果你不使用浮點(diǎn)數(shù)，答案會(huì)更好（主觀地）。此外，我發(fā)現(xiàn)將變量保留x為函數(shù)x而不是應(yīng)用函數(shù)更自然x(t)，例如：
In [15]: x = Function('x')
In [16]: x
Out[16]: x
In [17]: x(t)
Out[17]: x(t)
In [18]: eq = x(t).diff(t, 2) + x(t).diff(t)/4 + 5*x(t)/2
In [19]: eq
Out[19]:
d
──(x(t)) 2
5?x(t) dt d
────── + ──────── + ───(x(t))
2 4 2
dt
In [20]: dsolve(eq, x(t), ics={x(1): 0, x(2): 5})
Out[20]:
? 1/4 ?√159?t? ? -t
?5?? ?sin?──────? ? ───
? ? 8 ? 1/4 ?√159?t?? 8
x(t) = ?────────────────── - 5?? ?cos?──────????
? ?√159? ? 8 ??
? tan?────? ?
? ? 8 ? ?
添加回答
舉報(bào)
0/150
提交
取消