2017-08-09 57 views
0

让我们用sympy找到线性系统的外力不能代替积分常数到响应的积分

from sympy import * 

t, w, beta = symbols('t omega beta', positive=1) 
x0, v0 = symbols('x0 v0') 
x = symbols('x', cls=Function) 
homogeneous = diff(x(t), t, 2)/w**2+x(t) 
force = sin(beta*w*t) 
disp = dsolve(homogeneous-force, x(t)).rhs 
display(disp) 
constants = solve((disp.subs(t,0)-x0, disp.diff(t).subs(t,0)-v0)) 
display(constants) 

获得

        sin(beta*omega*t) 
C1*sin(omega*t) + C2*cos(omega*t) - ----------------- 
              2   
             beta - 1  
      2        
     beta *v0 + beta*omega - v0   
[{C1: --------------------------, C2: x0}] 
       / 2 \    
      omega*\beta - 1/    

现在,来完成我的解决方案的响应,我试图用这些常数的积分代入一般积分+积分

display(disp.subs(constants)) 

,给我

        sin(beta*omega*t) 
C2*sin(omega*t) + C2*cos(omega*t) - ----------------- 
              2   
             beta - 1  

我在哪里犯过错误?

+0

什么是你期望的输出? – jmoon

+1

@jmoon - 在第二块代码中,我希望看到代替“C1”和“C2”的符号值,这些值是我在第一块代码块的末尾计算和显示的,因为我已经让sympy去操作这样的替换。 – gboffi

回答

2

subs的参数应该是一个字典。

你从solve得到的不是字典,而是包含字典的单元素列表。使用它作为

disp.subs(constants[0]) 

那么结果是

x0*cos(omega*t) - sin(beta*omega*t)/(beta**2 - 1) + (beta**2*v0 + beta*omega - v0)*sin(omega*t)/(omega*(beta**2 - 1)) 
+0

哎呀,我知道它更好......我不得不说我被我之前发现的[答案](https://stackoverflow.com/a/41193705/2749397)部分误导了我。非常感谢这两个答案。 – gboffi