常微分方程组求数值解python+matlab

python实现

常微分函数表达式：

import scipy.integrate as spi
import numpy as np
from matplotlib import pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# 参数
ATR=0.2
ATAR=0.3
SR=0.1
SAR=0.3
SRC=0.7
PRC=0.2
PNP=0.3
AN=0
VR=SR+ATR+AN
VAR=SAR+ATAR+PRC*SRC
VI=PNP
a12=(VR-VAR)/3
a21=(-VR+VAR)/3
a31=(-VR+VI)/3
a32=(-VAR+VI)/3
RS=0.25
ARS=0.1
I=0.65
y=[RS,ARS,I]
k=9
w1=(k*k+k)/((k+3)*(k-2))#2
w2=3/((k+3)*(k-2))#0.5
w3=(k*k+k-3)/((k+3)*(k-2))
def diff_eqs(y,t,w1,w2,w3,a12,a31,a32,a21):
    V=y
    X=np.zeros((3))
    X[0] = V[0]*(w1*V[1]*a12-w2*V[2]*a31-V[0]*V[2]*a31-V[1]*V[2]*a32)
    X[1] = V[1]*(w1*V[0]*a21-w2*V[2]*a32-V[0]*V[2]*a31-V[1]*V[2]*a32)
    X[2] = V[2]*(w3*V[0]*a31+w3*V[1]*a32-V[0]*V[2]*a31-V[1]*V[2]*a32)
    return X
t_range=np.arange(0,200,0.1)
y1=[0.25,0.1,0.65]
# 求常微分在每个t时刻的值
RES= spi.odeint(diff_eqs,y1,t_range,args=(w1,w2,w3,a12,a31,a32,a21))
plt.plot(t_range,RES[:,0], '-r', marker='*', label='RS')
plt.plot(t_range,RES[:,1], '-g', marker='o', label='ARS')#换另一个可能是增加的
plt.plot(t_range,RES[:,2], '-b', marker='x', label='I')
plt.legend(loc=0)
plt.xlabel('iteration')
plt.ylabel('RS/ARS/I fraction')
plt.ylim(0,1)
plt.xlim(0,200)
plt.show()

matlab实现

function dydt= odefun( t,y )
ATR=0.2;
ATAR=0.3;
SR=0.1;
SAR=0.3;
SRC=0.7;
PRC=0.2;
PNP=0.3;
AN=0;
VR=SR+ATR+AN;
VAR=SAR+ATAR+PRC*SRC;
VI=PNP;
k=3;
w1=(k*k+k)/((k+3)*(k-2));
w2=3/((k+3)*(k-2));
w3=(k*k+k-3)/((k+3)*(k-2));
a12=(VR-VAR)/3;
a21=(-VR+VAR)/3;
a31=(-VR+VI)/3;
a32=(-VAR+VI)/3;
dydt=zeros(3,1);
dydt(1)=y(1)*(w1*y(2)*a12-w2*y(3)*a31-y(1)*y(3)*a31-y(2)*y(3)*a32);
dydt(2)=y(2)*(w1*y(1)*a21-w2*y(3)*a32-y(1)*y(3)*a31-y(2)*y(3)*a32);
dydt(3)=y(3)*(w3*y(1)*a31+w3*y(2)*a32-y(1)*y(3)*a31-y(2)*y(3)*a32);
end

函数调用：

tspan=0:1:200;
[t,y]=ode45('odefun',[0,200],[0.25 0.1 0.65]);
plot(t,y(:,1),'-o',t,y(:,2),t,y(:,3))
set(gca,'YLim',[0,1])
set(gca,'XLim',[0,200])