Index in position 2 exceeds array bounds (Index must not exceed 1) errors?

4 views (last 30 days)
JEESHAN ANSARI on 16 Jan 2024
Commented: JEESHAN ANSARI on 19 Jan 2024
% Code to solve system of four simultaneous differential equation
clc
N0 = [9.66e17;0;0;0]; % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
[T, N] = ode15s(@rate_eq,[0 126e-15], N0); % Intital conditions of time from 0 to 126e-15
plot(T,N(:,4),T,N(:,3),T,N(:,2),T,N(:,1)) % Plotting of N1 N2 N3 N4 vs Time
function dN = rate_eq( t,N )
t = 126e-15;
W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;
dN = zeros(4,7); % Blank array to store N1 N2 N3 N4 from solution for each values of seven values of W and PION
Q = 10e10; % Constant value
A = 10e8; % Constant value
Ppd = 10e10; % Constant value
% W and PION is not constant and its values are calculated in first part of
% code using FOR loop
for j = 1:length(W)
dN(1,j) = (-W(j)*N(1,j)) + ((Q+A)*N(2,j)) + ((W(j)+Q+A)*N(3,j)) + 0; % first equation
dN(2,j) = 0 + (-(Q+A)*N(2,j)) + ((Q+A)*N(3,j)) + 0; % second equation
dN(3,j) = (W(j)*N(1,j)) + 0 - (W(j)+(2*Q)+(2*A)+ Ppd + PION(j))*N(3,j) + 0; % third equation
dN(4,j) = 0 + 0 + PION(j)*N(3,j) + 0; % fourth equation
end
end
%-----------end of code----------------
I want to store values of N1 N2 N3 N4 for each seven different values of W and PION.
I also want to plot N1 N2 N3 N4 vs time for each W and PION.
But I am getting the following errors:
Index in position 2 exceeds array bounds. Index must not exceed 1.
Error in test_soln_rate>rate_eq (line 24)
dN(1,j) = (-W(j)*N(1,j)) + ((Q+A)*N(2,j)) + ((W(j)+Q+A)*N(3,j)) + 0; % first equation
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 148)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in test_soln_rate (line 6)
[T, N] = ode15s(@rate_eq,[0 126e-15], N0); % Intital conditions of time from 0 to 126e-15

Alan Stevens on 16 Jan 2024
N0 is a 4x1 vector, but you call it with two indices in the for loop! Turn N(1,j) into N(1) and similarly for the other calls to N.
Alan Stevens on 19 Jan 2024
Like this?
N0 = [9.66e17;0;0;0]; % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
dt = 126e-15/50;
tspan = 0:dt:126e-15;
W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;
for j = 1:7
w = W(j); p = PION(j);
[T, N] = ode15s(@(t,N) rate_eq(t,N,w,p),tspan, N0);
N1(:,j) = N(:,1); N2(:,j) = N(:,2); N3(:,j) = N(:,3); N4(:,j) = N(:,4);
end
% Plotting of N1 N2 N3 N4 vs Time
figure
plot(T,N1),grid
xlabel('T'),ylabel('N1')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','southwest')
figure
plot(T,N2),grid
xlabel('T'),ylabel('N2')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')
figure
plot(T,N3),grid
xlabel('T'),ylabel('N3')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')
figure
plot(T,N4),grid
xlabel('T'),ylabel('N4')
legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7','Location','northwest')
function dN = rate_eq(~,N,w,p )
Q = 10e10; % Constant value
A = 10e8; % Constant value
Ppd = 10e10; % Constant value
N1 = N(1); N2 = N(2); N3 = N(3); N4 = N(4);
dN1 = (-w*N1) + ((Q+A)*N2) + ((w+Q+A)*N3) + 0; % first equation
dN2 = 0 + (-(Q+A)*N2) + ((Q+A)*N3) + 0; % second equation
dN3 = (w*N1) + 0 - (w+(2*Q)+(2*A)+ Ppd + p)*N3 + 0; % third equation
dN4 = 0 + 0 + p*N3 + 0; % fourth equation
dN = [dN1; dN2; dN3; dN4];
end
JEESHAN ANSARI on 19 Jan 2024
Thanks @Alan Stevens for updating the code for me again. Now I am getting N1 N2 N3 N4 arrays as desired.

Categories

Find more on Data Type Identification in Help Center and File Exchange

R2023b

Community Treasure Hunt

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

Start Hunting!