1.问题：求系统的零状态响应

image.png

2.引入

首先用高数知识求解非齐次常系数微分方程

image.png

再利用信号与系统中冲激响应求解验证

image.png

利用MATLAB求解验证

y=dsolve('D2y+3*Dy+2*y=exp(-t)','y(0)=1','Dy(0)=2','t')

得出结果：

y =

(t - 2 exp(-t) + 3) exp(-t)

根据结果检验，上述手动计算与实际计算机得出结果一致。

t=0:0.1:20;

y = (t - 2 .*exp(-t) + 3) .*exp(-t);

y1=-exp(-t) .*(t - 2 .*exp(-t) + 3) + exp(-t).* (1 + 2.* exp(-t));

plot(t,y,'r-',t,y1,'b-'),legend('y','y’')

用MATLAB模拟图像结果：

image.png

3.利用Python求解该方程

通过上述计算，我们利用Python求解系统的零状态响应：

库函数准备

scipy

sympy

matplotlib

numpy

利用sympy进行符号解法

from sympy import *

init_printing()

#定义符号常量x 与 f(x)

x = Symbol('x')

f = symbols('f', cls=Function)

#用diffeq代表微分方程： f''(x) + 3f'(x) + 2f(x) = exp(-x)

diffeq = Eq(f(x).diff(x, x) + 3*f(x).diff(x) + 2*f(x), exp(-x))

#调用dsolve函数,返回一个Eq对象，hint控制精度

print(dsolve(diffeq, f(x)))

得到符号解，输出如下

Eq(f(x), (C1 + C2*exp(-x) + x)*exp(-x))

在带入初始松弛条件：

C1=-2

C2=3

结果与我们计算结果一致。

利用Numpy和Scipy进行数值解法

import matplotlib.pyplot as plt

from scipy import linspace,exp

from scipy.integrate import odeint, solve_bvp, solve_ivp

import numpy as np

'''

为了兼容solve_ivp的参数形式，微分方程函数定义的参数顺序为(t,y)，因此使用odeint函数时需要使参数tfirst=True

二阶甚至高阶微分方程组都可以变量替换成一阶方程组的形式，再调用相关函数进行求解，因此编写函数的时候，不同于一阶微分方程，二阶或者高阶微分方程返回的是低阶到高阶组成的方程组，

'''

def fvdp1(t,y):

'''

要把y看出一个向量，y = [dy0,dy1,dy2,...]分别表示y的n阶导，那么

y[0]就是需要求解的函数，y[1]表示一阶导，y[2]表示二阶导，以此类推

'''

dy1 = y[1] # y[1]=dy/dt，一阶导

dy2 = -3 * y[1] - 2 * y[0] + exp( -1 * t )

# y[0]是最初始，也就是需要求解的函数

# 注意返回的顺序是[一阶导， 二阶导]，这就形成了一阶微分方程组

return [dy1,dy2]

# 或者下面写法更加简单

def fvdp2(t,y):

'''

要把y看出一个向量，y = [dy0,dy1,dy2,...]分别表示y的n阶导

对于二阶微分方程，肯定是由0阶和1阶函数组合而成的，所以下面把y看成向量的话，y0表示最初始的函数，也就是我们要求解的函数，y1表示一阶导，对于高阶微分方程也可以以此类推

'''

y0, y1 = y

# y0是需要求解的函数，y1是一阶导

# 返回的顺序是[一阶导， 二阶导]，这就形成了一阶微分方程组

dydt = [y1, -3*y1-2*y0+exp(-t)]

return dydt

def solve_second_order_ode():

'''

求解二阶ODE

'''

t2 = linspace(0,20,1000)

tspan = (0, 20.0)

y0 = [1.0, 2.0] # 初值条件

# 初值[2,0]表示y(0)=2,y'(0)=0

# 返回y，其中y[:,0]是y[0]的值，就是最终解，y[:,1]是y'(x)的值

y = odeint(fvdp1, y0, t2, tfirst=True)

y_ = solve_ivp(fvdp2, t_span=tspan, y0=y0, t_eval=t2)

plt.subplot(211)

y1, = plt.plot(t2,y[:,0],label='y')

y1_1, = plt.plot(t2,y[:,1],label='y‘')

plt.legend(handles=[y1,y1_1])

plt.subplot(212)

y2, = plt.plot(y_.t, y_.y[0,:],'g--',label='y(0)')

y2_2, = plt.plot(y_.t, y_.y[1,:],'r-',label='y(1)')

plt.legend(handles=[y2,y2_2])

plt.show()

solve_second_order_ode()

结果：

image.png