sir.m

这是SIR微分方程模型对初值问题的求解和作图

global sigma Y0  Y1 options; %全局变量
Y0=[0.02,0.7]; Y1=[0.05,0.6]; 
sigma=log(Y1(2)/Y0(2))/((Y1(1)+Y1(2))-(Y0(1)+Y0(2)))
opt=optimset('TolX',1e-10,'TolFun',1e-20); %自设定以下误差
[mu, fnorm]=lsqnonlin(@fun, 0.3,[],[],opt);
disp(['mu = ',num2str(mu,'%12.8f'), ', ', 'sigma = ',num2str(sigma,'%12.8f')]);
options=odeset( 'RelTol',4e-14,'AbsTol', [1e-15,1e-15]); %自定义常微分方程的解的精度
[T,Y]=ode45(@odefun, [0,50], Y0, options, mu);
disp(['i(50)=', num2str(Y(end,1),'%12.8f'),',  ',' s(50)=', num2str(Y(end,2),'%12.8f')]);
figure; %打开新的图形窗口
plot(Y(:,2),Y(:,1)); %画相轨线
title('s-i 相轨线'); %写标题
xlabel('s(t)'); ylabel('i(t)'); %标轴的名称
figure; %打开新的图形窗口
plot(T, Y(:,2),T,Y(:,1));
title('s-i 图像'); %写标题
xlabel('t'); ylabel('i(t),s(t)'); %标轴的名称
text(T(100), Y(100,1), 'i(t)','FontSize', 16); %标注
text(T(100), Y(100,2), 's(t)','FontSize', 16); %标注
legend('s(t)','i(t)'); %图例