美赛 4：微分方程（编程篇）

Destiny坠明

已于 2022-01-26 01:53:59 修改

阅读量2.5k

点赞数

分类专栏：北师大数学模型文章标签： matlab python

于 2022-01-26 01:49:18 首次发布

本文链接：https://blog.csdn.net/qq_55818063/article/details/122682259

版权

北师大数学模型专栏收录该内容

15 篇文章 170 订阅

订阅专栏

一、Matlab & Python 求解微分方程

1.求方程解析解

2.求方程数值解（绘图）

二、颜色盘

一、Matlab & Python 求解微分方程

1.求方程解析解

例1：dy/dx = 2x , y(1) = 2

Matlab解法：

dsolve('y-Dy=2*x','y(0)=3','x')

Python解法：

import sympy as sy

x=sy.symbols('x')
y=sy.symbols('y',cls=sy.Function)

diffeq=sy.Eq(y(x).diff(x),2*x)

print(sy.dsolve(diffeq,y(x),ics={y(1):2}))
#sy.pprint(sy.dsolve(diffeq,y(x),ics={y(1):2}))

例2：d2s/dt2 = -0.4 , s(0) = 0 , Ds(0) = 20

Matlab解法：

dsolve('D2s=-0.4','s(0)=0,Ds(0)=20','t')

Python解法：

t=sy.symbols('t')
s=sy.symbols('s',cls=sy.Function)
a=sy.symbols('-0.4')

diffeq=sy.Eq(s(t).diff(t,2),a)
print(sy.dsolve(diffeq,s(t),ics={s(0):0,s(t).diff(t).subs(t,0):20}))

例3：dy/dx = 2xy

Matlab解法：

dsolve('Dy=2*x*y','x')

Python解法：

x=sy.symbols('x')
y=sy.symbols('y',cls=sy.Function)

diffeq=sy.Eq(y(x).diff(x),2*x*y(x))
print(sy.dsolve(diffeq,y(x)))

例4：dM/dt = -nM , M(0) = M0

Matlab解法：

dsolve('DM=-n*M','M(0)=M0','t')

Python解法：

t=sy.symbols('t')
n=sy.symbols('n')
M0=sy.symbols('M0')
M=sy.symbols('M',cls=sy.Function)
#Y=sy.symbols('Y',cls=sy.Function)
#a=sy.symbols('-0.4')

diffeq=sy.Eq(M(t).diff(t),-n*M(t))
print(sy.dsolve(diffeq,M(t),ics={M(0):M0}))

例5：m*dv/dt = mg-kv , v(0) = 0

Matlab解法：

dsolve('m*Dv=m*g-k*v','v(0)=0','t')

Python解法：

t=sy.symbols('t')
m=sy.symbols('m')
g=sy.symbols('g')
k=sy.symbols('k')
v=sy.symbols('v',cls=sy.Function)
#Y=sy.symbols('Y',cls=sy.Function)
#a=sy.symbols('-0.4')

diffeq=sy.Eq(m*v(t).diff(t),m*g-k*v(t))
print(sy.dsolve(diffeq,v(t),ics={v(0):0}))

例6：dt/dh*KS（2gh）**(1/2) = -pi(2h-h**2) , t(1) = 0

Matlab解法：

dsolve('Dt*K*S*(2*g*h)^(1/2)=-pi*(2*h-h^2)','t(1)=0','h')

Python解法：

pi=sy.symbols('pi')
g=sy.symbols('g')
h=sy.symbols('h')
S=sy.symbols('S')
K=sy.symbols('K')
t=sy.symbols('t',cls=sy.Function)

diffeq=sy.Eq(t(h).diff(h)*K*S*(2*g*h)**0.5,-pi*(2*h-h**2))
print(sy.dsolve(diffeq,t(h),ics={t(1):0}))

快速定义并计算数学函数：

f = lambda x : function(x)

例7：y**2+x**2*dy/dx = xy*dy/dx

Matlab解法：

dsolve('y^2+x^2*Dy=x*y*Dy','x')

Python解法：

x=sy.symbols('x')
y=sy.symbols('y',cls=sy.Function)

diffeq=sy.Eq(y(x)**2+y(x).diff(x)*x**2,x*y(x)*y(x).diff(x))
print(sy.dsolve(diffeq,y(x)))

例8：y/y' -x = (x**2+y**2)**(1/2)

Matlab解法：

dsolve('Dy-2*y/(x+1)=(x+1)^(5/2)','x')

Python解法：

x=sy.symbols('x')
y=sy.symbols('y',cls=sy.Function)

diffeq=sy.Eq(y(x).diff(x)-2*y(x)/(x+1),(x+1)**(5/2))
print(sy.dsolve(diffeq,y(x),ics={y(0):1}))

y1 = lambda x:(x+1)**2 * (2/3 * (x+1)**(3/2) + 1/3)
y2 =lambda x:(0.333333333333333*x**4 + 1.33333333333333*x**3 + 2.0*x**2 + 1.33333333333333*x + 0.666666666666667*(x + 1)**5.5 + 0.333333333333333)/(1.0*x**2 + 2.0*x + 1.0)

# y1 == y2

例9：y''' = e**2x - cos(x)

Matlab解法：

dsolve('D3y=exp(2*x)-cos(x)','x')

Python解法：

x=sy.symbols('x')
y=sy.symbols('y',cls=sy.Function)

diffeq=sy.Eq(y(x).diff(x,3),sy.exp(2*x)-sy.cos(x))
print(sy.dsolve(diffeq,y(x)))

总结：求解析解的问题，Matlab代码简洁，运算效率高；

Python仅作参考。。。

2.求方程数值解（绘图）

例：改进的SEIR传染病模型

Matlab解法：

%% SEIR.m

function dy=SEIR(t,y)
    
    u=0.002;
    beta1=0.1;beta2=0.05;
    sigma=0.2;
    gamma1=0.01;gamma2=0.01;
    v=0.001;
    d=0.005;
    alpha=0.01;

    dy=zeros(7,1);
    S=y(1);
    E=y(2);
    I=y(3);
    R1=y(4);
    R2=y(5);
    ID=y(6);
    ND=y(7);

    N=y(1)+y(2)+y(3)+y(4)+y(5);

    dy(1)=-beta1*S*I/N-beta2*S*E/N+alpha*R2+u*N-v*S;
    dy(2)=beta1*S*I/N+beta2*S*E/N-sigma*E-v*E;
    dy(3)=sigma*E-d*I-gamma1*I-gamma2*I-v*I;
    dy(4)=gamma1*I-v*R1;
    dy(5)=gamma2*I-v*R2-alpha*R2;
    dy(6)=d*I;
    dy(7)=v*N;

end

%% main.m

[t,y]=ode45('SEIR',[0:500],[1000,1,1,0,0,0,0]);
hold on;
plot(t,y(:,1),'LineWidth',1.5,'Color',[240/255,92/255,30/255]);
plot(t,y(:,2),'Color',[255/255,186/255,17/255]);
plot(t,y(:,3),'Color',[235/255,24/255,56/255]);
plot(t,y(:,4),'Color',[4/255,160/255,122/255]);
plot(t,y(:,5),'Color',[123/255,188/255,70/255]);
plot(t,y(:,6),'Color',[159/255,34/255,134/255]);
plot(t,y(:,7),'Color',[8/255,76/255,161/255]);

legend('S','E','I','R1','R2','ID','ND');

Python解法：

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

def SEIRS(y,t,u,beta1,beta2,sigma,gamma1,gamma2,v,d,alpha):
    S,E,I,R1,R2,ID,ND = y
    N=S+E+I+R1+R2

    return np.array([-beta1*S*I/N-beta2*S*E/N+alpha*R2+u*N-v*S,
                     beta1*S*I/N+beta2*S*E/N-sigma*E-v*E,
                     sigma*E-d*I-gamma1*I-gamma2*I-v*I,
                     gamma1*I-v*R1,
                     gamma2*I-v*R2-alpha*R2,
                     d*I,
                     v*N])

t=np.linspace(0,500,50000)

y=odeint(SEIRS,[1000,1,1,0,0,0,0],t,args=(0.002,0.1,0.05,0.2,0.01,0.01,0.001,0.005,0.01))

I_max_y=np.max(y[:,2])
I_max_t=t[np.argmax(y[:,2])]

plt.plot(t,y[:,0],label='S',color=(240/255,92/255,30/255))
plt.plot(t,y[:,1],label='E',color=(255/255,186/255,17/255))
plt.plot(t,y[:,2],label='I',lw=2.5,color=(235/255,24/255,56/255))
plt.plot(t,y[:,3],label='R1',color=(4/255,160/255,122/255))
plt.plot(t,y[:,4],label='R2',color=(123/255,188/255,70/255))
#plt.plot(t,y[:,5],label='ID',color=(159/255,34/255,134/255))
#plt.plot(t,y[:,6],label='ND',color=(8/255,76/255,161/255))

plt.grid(linestyle=":")
plt.legend()

plt.plot(I_max_t,I_max_y,'.',markersize=10,color=(235/255,24/255,56/255))

plt.annotate("maximum",xy=(I_max_t+1,I_max_y+5),xytext=(172,580),weight="bold",color=(235/255,24/255,56/255),
             arrowprops=dict(arrowstyle="fancy",connectionstyle="arc3",color=(235/255,24/255,56/255)))

#plt.annotate("",xy=(I_max_t-8.5,I_max_y+5),xytext=(12,10),color=(235/255,24/255,56/255),
#             arrowprops=dict(arrowstyle="simple",color=(235/255,24/255,56/255)))

#plt.annotate("",xy=(266,135),xytext=(I_max_t+13.5,I_max_y+3.5),color=(235/255,24/255,56/255),
#             arrowprops=dict(arrowstyle="simple",connectionstyle="arc3",color=(235/255,24/255,56/255)))

plt.show()