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

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

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

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

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

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def f(z,t,u):
    dzdt = np.zeros(2,)
    x = z[0]
    y = z[1]
    dzdt[0] = (-x + u) / 2.0
    dzdt[1] = (-y + x) / 5.0
    return dzdt

z0  = [0,0]
n = 150
t = np.linspace(0,15, n)
u = np.zeros(n)
u[0] = 1
x = np.zeros(n)
y = np.zeros(n)
for i in range (1,n):
    tspan = [t[i-1],t[i]]
    z = odeint(f,z0,t,args=(u[i],))
    u[i]  = x[i] + 3*y[i]
    z0 = z[1]
    x[i] = z0[0]
    y[i] = z0[1]

plt.plot(t,u)
plt.plot(t,x)
plt.plot(t,y)
plt.show()
kiayqfof

kiayqfof1#

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

In [34]: x, y, u = symbols('x, y, u', cls=Function)

In [35]: t = symbols('t')

In [36]: eqs = [Eq(x(t).diff(t), (-x(t) + u(t)) / 2.0),
    ...:        Eq(y(t).diff(t), (-y(t) + x(t)) / 5.0),
    ...:        Eq(u(t), x(t) + 3*y(t))]

In [37]: eqs[2]
Out[37]: u(t) = x(t) + 3⋅y(t)

In [38]: eqs_new = [eq.subs(eqs[2].lhs, eqs[2].rhs) for eq in eqs[:2]]

In [39]: eqs_new
Out[39]: 
⎡d                    d                             ⎤
⎢──(x(t)) = 1.5⋅y(t), ──(y(t)) = 0.2⋅x(t) - 0.2⋅y(t)⎥
⎣dt                   dt                            ⎦

In [41]: dsolve(eqs_new)[0]
Out[41]: 
                              -0.656776436283002⋅t                        0.456776436283002⋅t
x(t) = - 2.28388218141501⋅C₁⋅ℯ                     + 3.28388218141501⋅C₂⋅ℯ                   

In [42]: dsolve(eqs_new)[1]
Out[42]: 
               -0.656776436283002⋅t           0.456776436283002⋅t
y(t) = 1.0⋅C₁⋅ℯ                     + 1.0⋅C₂⋅ℯ

相关问题