[未答复] matlab求解微分方程画图改变参数取值图就出不来了,求解!

[复制链接]
violet1024 发表于 2022-1-12 10:05:36
function dy=shi(~,y,~)
A=15;
beta=0.02;
beta_w=0.2143;
gamma=0.004;
alpha=10;
theta=0.003;
phi=0.4;
xi=0.33;
varepsilon=0.2;
mu=0.0000548;
N=80;
a1=1/(gamma+mu);
a2=alpha/(a1*xi);
sum_k = 0;
for k=1:N
    r_k=(exp(-40)*(40^k)/factorial(k));
    N_k=A*r_k/mu;
    s_k=(theta+mu)*N_k/(phi+theta+mu);
    v_k=phi*N_k/(phi+theta+mu);
    sum_k = sum_k + k*N_k;
end
    sum_k=sum(sum_k);
C=zeros(80,80);
for k=1:N
for j=1:N
    tau_k=beta*k*s_k/sum_k;
    delta_k=varepsilon*beta*k*v_k/sum_k;
    sigma_k=beta_w*s_k;
C(k,j)=j*a1*(tau_k+delta_k)+a2*sigma_k;
end
end
R_0=vrho(C)
l = zeros(1,80);
for i=1:80
    l(i) = i;
end
kP=l;
dy=zeros(3*N+1,1);
    sum_m=kP*(y(160+1:240))
    for j=1:N
        
        dy(j)=A*r_k-beta*j*y(j)*(sum_m/sum_k)-beta_w*y(j)*y(241)+theta*y(80+j)-(phi+mu)*y(j);
        
        dy(80+j)=phi*y(j)-varepsilon*beta*j*y(j)*(sum_m/sum_k)-(theta+mu)*y(80+j);
        
        dy(160+j)=beta*j*y(j)*(sum_m/sum_k)+varepsilon*beta*j*y(j)*(sum_m/sum_k)+beta_w*y(j)*y(241)-(gamma+mu)*y(160+j);


        dy(241)=alpha*sum_m-xi*y(241);
    end
end




clear;
clc;
k_max=80;
    I1=zeros(1,k_max);
    I2=zeros(1,k_max);
    I3=zeros(1,k_max);
    I4=zeros(1,1);
for j=1:80
    I1(j)=500;
    I2(j)=100;
    I3(j)=100;
end
I4(1)=300;
y0=[I1';I2';I3';I4'];
% y1 = y0';
disp('画图前');
tspan = (0:1:150);
[T,Y]=ode45(@shi,tspan,y0);
disp('画图');
plot(T,Y(:,200),'b','linewidth',1.5)
只是改变最上面的参数取值,图就出不来了,陷入了死循环,求大神帮忙看下哪里有问题,感谢!

您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表