1 view (last 30 days)

Show older comments

I have been working on this for a while now, but can't seem to figure out what's the problem. I keep getting the error:

Attempted to access betaI(2.72775e+06); index out of bounds because numel(betaI)=1.

for the following code in MATLAB: (Function file followed by script)

%System of differential equations

function ydot = D12Eqns(t,y)

global Lambda betaI betaA tau mu deltaIL nuIL muA deltaAL muT r2 r1 nuIS deltaSL nuIT h1 h2 a1 a2

XSS = y(1);

XSL = y(2);

XST = y(3);

XIS = y(4);

XIL = y(5);

XIT = y(6);

XAS = y(7);

XAL = y(8);

XAT = y(9);

TSC = y(10);

TIC = y(11);

TAC = y(12);

N = XSS+XSL+XST+XIS+XIL+XIT+XAS+XAL+XAT+TSC+TIC+TAC;

ydot = zeros(12,1);

ydot(1) = Lambda-betaI(XIS+XIL+XIT)*XSS/N-betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XSS/N-mu*XSS;

ydot(2) = tau(XST+XIT+XAT)*XSS/N-betaI(XIS+XIL+XIT)*XSL/N-betaA(XAS+XAL+XAT)*XSL/N-(deltaSL+mu)*XSL-r1*XSL;

ydot(3) = deltaSL*XSL-betaI(XIS+XIL+XIT)*XST/N-betaA(XAS+XAL+XAT)*XST/N-(mu+mu*T)*XST-r2*XST;

ydot(4) = betaI(XIS+XIL+XIT)*XSS/N+betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XIS/N-(nuIS+mu)*XIS;

ydot(5) = betaI(XIS+XIL+XIT)*XSL/N+betaA(XAS+XAL+XAT)*XSL/N+tau(XST+XIT+XAT)*XIS/N-(deltaIL+nuIL+mu)*XIL-h1*XSL;

ydot(6) = betaI(XIS+XIL+XIT)*XST/N+betaA(XAS+XAL+XAT)*XST/N+deltaIL*XIL-(nuIT+mu+muT)*XIT-h2*XIT;

ydot(7) = nuIS*XIS-tau(XST+XIT+XAT)*XAS/N-(mu+muA)*XAS;

ydot(8) = nuIL*XIL+tau(XST+XIT+XAT)*XAS/N-(deltaAL+mu+muA)*XAL-a1*XAL;

ydot(9) = nuIT*XIT+deltaAL*XAL-(mu+muA+muT)*XAT-a2*XAT;

ydot(10)= r1*XSL+r2*XST-mu*TSC-betaI(XIS+XIL+XIT)*TSC/N-betaA(XAS+XAL+XAT)*TSC/N;

ydot(11)= h1*XSL+h2*XIT-nuIS*TIC-mu*TIC;

ydot(12)= a1*XAL+nuIS*TIC+a2*XAT-(muA+muT)*TAC;

%Project Model Solver

clear all

global Lambda betaI betaA tau mu deltaIL nuIL muA deltaAL muT r2 r1 nuIS deltaSL nuIT h1 h2 a1 a2

Lambda = 6336000; %0.022 x N = 0.022 x 288000000%

betaI = 0.170;

betaA = 0.204;

tau = 4;

mu = 0.016;

deltaIL = 0.05;

muA = 0.5;

muT = 0.2;

deltaAL = 0.1;

deltaSL = 0.0038;

nuIS = 0.085;

nuIT = 1;

nuIL = 1.7;

r1 = 3;

r2 = 1;

h1 = 3;

h2 = 12;

a1 = 12;

a2 = 12;

%%%== Initial conditions ==

XSS0 = 168558250;

XSL0 = 112668000; %XSL0-50000

XST0 = 1968750; %XST0-50000

XIS0 = 1587563;

XIL0 = 1085500; %XIL0-25000

XIT0 = 54687; %XIT0-25000

XAS0 = 529188; %XAS0-12500

XAL0 = 358000; %XAL0-12500

XAT0 = 14062;

TSC0 = 100000; %In TSC0, TIC0, and TAC0, we should reduce the numbers in the rest of the compartments (lines 27 to 35) and subtract from the total population and then divide them among the three compartments.

TIC0 = 50000;

TAC0 = 25000;

y0 = [XSS0 XSL0 XST0 XIS0 XIL0 XIT0 XAS0 XAL0 XAT0 TSC0 TIC0 TAC0];

tf = 40; %years;

[t, y] = ode15s(@D12Eqns,[0 tf],y0);

XSS = y(:,1);

XSL = y(:,2);

XST = y(:,3);

XIS = y(:,4);

XIL = y(:,5);

XIT = y(:,6);

XAS = y(:,7);

XAL = y(:,8);

XAT = y(:,9);

TSC = y(:,10);

TIC = y(:,11);

TAC = y(:,12);

plot(t,100.*(TSC)./N,'b')

hold on

axis square

xlabel('Years')

ylabel('Percentage of people in TSC out of whole population')

Would someone mind helping me? I'm working with a system of 12 ordinary differential equations. Thank you for your time and patience.

Laurent
on 31 Oct 2013

Without running your code, I see that you define betaI as a scalar:

betaI = 0.170;

But then in your function you use betaI as if it is an array/vector:

ydot(1) = Lambda-betaI(XIS+XIL+XIT)*XSS/N-betaA(XAS+XAL+XAT)*XSS/N-tau(XST+XIT+XAT)*XSS/N-mu*XSS;

and so on.

You are trying to use the 'XIS+XIL+XIT'th value of betaI, while there is only 1, namely 0.170.

So either betaI is not correctly defined, or your ydot(1) and all the others that use betaI are not calculated correctly. Did you maybe forget a * between betaI and (XIS+XIL+XIT)?

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!