scipy 用Python求解微分方程组和线性方程组

kmynzznz  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(241)

我试着解这个方程组,它是两个常微分方程和一个线性方程

  1. dx(t)/dt = (-x(t) + u(t)) / 2.0,
  2. dy(t)/dt = (-y(t) + x(t)) / 5.0
  3. u(t) = x(t) + 3*y(t)

在这个初始条件下,我一直在尝试用odeint来解这个方程,但是我不知道如何定义线性函数,也就是函数u,我不知道这是不是解这个方程的另一种方法。

  1. import numpy as np
  2. from scipy.integrate import odeint
  3. import matplotlib.pyplot as plt
  4. def f(z,t,u):
  5. dzdt = np.zeros(2,)
  6. x = z[0]
  7. y = z[1]
  8. dzdt[0] = (-x + u) / 2.0
  9. dzdt[1] = (-y + x) / 5.0
  10. return dzdt
  11. z0 = [0,0]
  12. n = 150
  13. t = np.linspace(0,15, n)
  14. u = np.zeros(n)
  15. u[0] = 1
  16. x = np.zeros(n)
  17. y = np.zeros(n)
  18. for i in range (1,n):
  19. tspan = [t[i-1],t[i]]
  20. z = odeint(f,z0,t,args=(u[i],))
  21. u[i] = x[i] + 3*y[i]
  22. z0 = z[1]
  23. x[i] = z0[0]
  24. y[i] = z0[1]
  25. plt.plot(t,u)
  26. plt.plot(t,x)
  27. plt.plot(t,y)
  28. plt.show()
kiayqfof

kiayqfof1#

你可以把u代入其他方程,这些特殊的常微分方程也可以用SymPy符号化求解:

  1. In [34]: x, y, u = symbols('x, y, u', cls=Function)
  2. In [35]: t = symbols('t')
  3. In [36]: eqs = [Eq(x(t).diff(t), (-x(t) + u(t)) / 2.0),
  4. ...: Eq(y(t).diff(t), (-y(t) + x(t)) / 5.0),
  5. ...: Eq(u(t), x(t) + 3*y(t))]
  6. In [37]: eqs[2]
  7. Out[37]: u(t) = x(t) + 3y(t)
  8. In [38]: eqs_new = [eq.subs(eqs[2].lhs, eqs[2].rhs) for eq in eqs[:2]]
  9. In [39]: eqs_new
  10. Out[39]:
  11. d d
  12. ⎢──(x(t)) = 1.5y(t), ──(y(t)) = 0.2x(t) - 0.2y(t)⎥
  13. dt dt
  14. In [41]: dsolve(eqs_new)[0]
  15. Out[41]:
  16. -0.656776436283002t 0.456776436283002t
  17. x(t) = - 2.28388218141501C₁⋅ℯ + 3.28388218141501C₂⋅ℯ
  18. In [42]: dsolve(eqs_new)[1]
  19. Out[42]:
  20. -0.656776436283002t 0.456776436283002t
  21. y(t) = 1.0C₁⋅ℯ + 1.0C₂⋅ℯ
展开查看全部

相关问题