目录
一、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 syx=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.mfunction 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 pltdef SEIRS(y,t,u,beta1,beta2,sigma,gamma1,gamma2,v,d,alpha):S,E,I,R1,R2,ID,ND = yN=S+E+I+R1+R2return 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()
总结：Matlab需要定义m文件来执行微分计算；
Python可以利用scipy库中的odeint进行数值求解，并结合matplotlib作图。
二、颜色盘
附：RGB颜色盘（QQ截图，按c复制色号） 注：归一化至 [0,1] 区间
