EDA365电子论坛网

标题: 索引超出矩阵维度怎么处理？ [打印本页]

作者: Touuqu 时间: 2020-8-12 15:20
标题: 索引超出矩阵维度怎么处理？
求解一个微分方程组，提示索引超出矩阵维度，但我看了好久也没找到哪里出问题了，请大神们帮帮忙。程序如下主程序：

clc
clear all;
close all;

tinit=0;                   % time min / initial time
tmax=100;                % time length
dt=0.01; % time interval
M=tmax/dt; % time points
tspan=linspace(tinit,M*dt,M+1);

K=100; %space points
I=zeros(K/2,1);
R=zeros(K/2,1);
T=zeros(1,1);
V=zeros(1,1);
%S=zeros(1,1);
%I=zeros(1,1);

s = 13;
d = 1;
beta =0.8;
c=20;
es=0;%without treatment es=0
ea=0;%without treatment ea=0
kappa=1;%%without treatment kappa=1

I(1)=beta;
i0=I;
R(K/2+1)=1; % gamma*I
r0=R;
%S(1)=S;
t0=1;
%I(1)=I;
v0=1;

aiinit=0.0;
aimax=84; % space length
dai=(aimax-aiinit)/K; %space interval
aispan=linspace(aiinit,aimax,K);
aispan=aispan.';

options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3],'AbsTol', 1e-4, 'NonNegative', [1 2 3]);
[t,sir] = ode15s(@(t,sir)mytest(t,sir,dai,aispan,K,beta,s,d,c,es,ea,kappa),tspan,[i0' r0' t0 v0],options);

figure
h1=surf(aispan,t,sir(:,1:K/2));
set(h1,'LineStyle','none');
hold on
xlabel('Age (ai)')
ylabel('Time (t)')
zlabel('I(ai,t)')

figure
h2=surf(aispan,t,sir(:,K/2+1:K));
set(h2,'LineStyle','none');
hold on
xlabel('Age (ai)')
ylabel('Time (t)')
zlabel('R(ai,t)')

figure
% R is matrix, sum(R,dimension) is a column vector containing sum of each row
plot(t,sum(sir(:,1:K/2),2),'Color',[0.5 0 0.9])
hold on
xlabel('Time (t)')
ylabel('I(t)')
legend('Recovered')

figure
% R is matrix, sum(R,dimension) is a column vector containing sum of each row
plot(t,sum(sir(:,K/2+1:K),2),'Color',[0.5 0 0.9])
hold on
xlabel('Time (t)')
ylabel('R(t)')
legend('Recovered')

figure
plot(t,sir(:,K+1),'Color',[0 0 1])
hold on
plot(t,sir(:,K+2))%,'Color',[1 0 0])
hold on
% R is matrix, sum(R,dimension) is a column vector containing sum of each row
plot(t,sum(sir(:,1:K/2),2),'Color',[0.5 0 0.9]);
hold on
plot(t,sum(sir(:,K/2+1:K),2),'Color',[0.5 0 0.9])
title('Population vs Time')
xlabel('Time (t)')
ylabel('Population')
legend('Susceptible','Infected','Recovered')

函数：
function dsir = mytest(t,sir,dai,aispan,K,beta,s,d,c,es,ea,kappa)

tau1=10;
tau2=15;
tau3=25;
tau4=35;

%rho=@(a) 0.*(a<=tau2)+(0.01*(a-tau2).^2.*exp(-0.2*(a-tau2))).*(a>tau2);
%delta=@(a) 0.*(a<=tau4)+(56^(-3)*(a-tau4).^2.*exp(-0.2*(a-tau4))).*(a>tau4);
%mu=@(a) 0.*(a<=tau3)+(22^(-3)*(a-tau3).^2.*exp(-0.2*(a-tau3))).*(a>tau3);
%alpha=@(a) 0.*(a<=tau1)+(0.2*(a-tau1).^2.*exp(-0.2*(a-tau1))).*(a>tau1);

ai=(0:K)*dai;
ai=ai ( : );
rho=zeros(numel(ai));
delta=zeros(numel(ai));
mu=zeros(numel(ai));
alpha=zeros(numel(ai));

I=sir(1:K/2);
R=sir(K/2+1:K);
T=sir(K+1);
V=sir(K+2);

dIda=zeros(K/2,1);
dRda=zeros(K/2,1);
for j=1:K
if (ai(j)<=tau2)
      rho(j)=0;
else
      rho(j)=0.01*(ai(j)-tau2)^2*exp(-0.2*(ai(j)-tau2));
end
if (ai(j)<=tau4)
      delta(j)=0;
else
      delta(j)=56^(-3)*(ai(j)-tau4)^2*exp(-0.2*(ai(j)-tau4));
end
if (ai(j)<=tau3)
      mu(j)=0;
else
      mu(j)=22^(-3)*(ai(j)-tau3)^2*exp(-0.2*(ai(j)-tau3));
end
if (ai(j)<=tau1)
      alpha(j)=0;
else
      alpha(j)=0.2*(ai(j)-tau1)^2*exp(-0.2*(ai(j)-tau1));
end
end

for j=1:K/2
if(j==1)
      dIda(j)=beta*V*T;
else
      dIda(j)=-1*(I(j)-I(j-1))/(aispan(j)-aispan(j-1))...
         -delta(j-1)*I(j-1);
end
end

for j=K/2+1:K
if(j==K/2+1)
      dRda(j)=1; % R(0,t)=gamma*I;
else
      dRda(j)=-1*(R(j)-R(j-1))/(aispan(j)-aispan(j-1))...
         +(1-ea)*alpha(j-1)-((1-es)*rho(j-1)+kappa*mu(j-1))*R(j-1);
end
end

ai_centered = (ai(1:numel(ai)-1)+ai(2:numel(ai)))/2.0;
value_integral = trapz(ai_centered,I.*R.*rho(ai_centered));
dTdt = s-d*T-beta*V*T;          % S=-beta*S*I;
dVdt = (1-es)*value_integral-c*V;    % I=beta*S*I-gamma*I;

dsir = [dIda;dRda;dTdt;dVdt];

end

一运行就报错：
索引超出矩阵维度。

出错 mytest (line 63)
      dRda(j)=-1*(R(j)-R(j-1))/(aispan(j)-aispan(j-1))...

出错 maintest>@(t,sir)mytest(t,sir,dai,aispan,K,beta,s,d,c,es,ea,kappa)

出错 odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

出错 ode15s (line 150)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

出错 maintest (line 43)
[t,sir] = ode15s(@(t,sir)mytest(t,sir,dai,aispan,K,beta,s,d,c,es,ea,kappa),tspan,[i0' r0' t0 v0],options);

作者: SsaaM7 时间: 2020-8-12 16:14
看看咯

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)