# I do not understand why my code still gives me errors

2 views (last 30 days)
DAVIDE TRONO on 23 Jul 2019
Edited: Jan on 24 Jul 2019
I'm basically setting a Runge-Kutta of order 4 to solve a sort of 'pred-prey' model. So, my three programs are:
1.
function [t,y]=RK(fun,tspan,y0,n,A,b,c)
h=(tspan(2)-tspan(1))/n;
y(:,1)=y0;
t(1)=tspan(1);
s=length(b);
for j=1:n
\$ k(:,1)=feval(fun,t(j),y(:,j));
for i=2:s
k(:,i)=feval(fun,t(j)+c(i)*h,y(:,j)+h*k(:,1:i-1)*A(i,1:i-1)');
end
end
2.
function dy=pred_prey(t,y)
function s=cos(x)
k=1;
a=2/3;
d=4/3;
r(t)=int(s(x^2),x,0,t);
mu(t)=13/20-(3/5)*exp(-(3/x));
dy(1)=(y(1)+k)*r(t)-a*y(1)*y(2);
dy(2)=-mu(t)*y(2)+d*y(1)*y(2);
dy=dy';
3.
function []=driver_pred_prey()
close all
tspan=[0 100];
n=1000;
fun='pred_prey';
y0=[0 -1]';
A=[0 0 0 0; 1/2 0 0 0 ; 0 1/2 0 0; 0 0 1 0];
b=[1/6 1/3 1/3 1/6]';
c=[0 1/2 1/2 1]';
\$\$ [t,y]=RK('pred_prey', tspan, y0, n, A, b, c);
figure (1)
plot(t, y(1,:), 'k-'), hold
plot(t, y(2,:), 'k-.')
legend('resources',ppl')
figure (2)
plot(y(1,:),y(2,:))
title('orbits)
but I really do not understand what errors there are at \$ and \$\$

DAVIDE TRONO on 23 Jul 2019
Well you are right, the errors reported by matlab are:
1. Output argument "dy" (and maybe others) not assigned during call to "pred_prey".
2. Error in RK (line 8) k(:,1)=feval(fun,t(j),y(:,j)); (this is the \$ line)
3. Error in driver_pred_prey (line 14) [t,y]=RK('pred_prey', tspan, y0, n, A, b, c);
the last one is the \$\$ line. To be clear, RK, pred_prey and driver_preded_prey are respectively the name of the codes 1,2 and 3.
I know that cos() is pretty easy to integrate, but cos(x^2) is not! Its primitive can't be written in terms of elementary functions.
Jan on 24 Jul 2019
% Integral cos2(x) dx = 0.5 * (cos(x)*sin(x) + x)
This is only the part of the messages, which shows where the error is, but not what it is:
• Error in RK (line 8) k(:,1)=feval(fun,t(j),y(:,j)); (this is the \$ line)
• Error in driver_pred_prey (line 14) [t,y]=RK('pred_prey', tspan, y0, n, A, b, c);
For the first message:
• Output argument "dy" (and maybe others) not assigned during call to "pred_prey".
Did you remove the orphaned "function s=cos(x)" as you have been instructed already? Then this message should disappear.
DAVIDE TRONO on 24 Jul 2019
function [t,y]=RK(fun,tspan,y0,n,A,b,c)
h=(tspan(2)-tspan(1))/n;
y(:,1)=y0;
t(1)=tspan(1);
s=length(b); %numero di l
for J=1:n
k(:,1)=feval(fun,t(j),y(:,j));
for i=2:s
k(:,i)=feval(fun,t(j)+c(i)*h,y(:,j)+h*k(:,1:i-1)*A(i,1:i-1)');
end
end
code number 1, it is called RK
function dy=pred_prey(t,y)
k=1;
a=2/3;
d=4/3;
f=@(x)cos(x.^2)
r=@(t)integral(f,0,t);
mu=@(t)13/20-(3/5)*exp(-(3/t));
dy(1)=(y(1)+k)*r(t)-a*y(1)*y(2);
dy(2)=-mu(t)*y(2)+d*y(1)*y(2);
dy=dy';
this is the number 2, called pred_prey
function []=driver_pred_prey()
close all
tspan=[0 500];
n=1000;
fun='pred_prey';
y0=[0 -1]';
A=[0 0 0 0; 1/2 0 0 0 ; 0 1/2 0 0; 0 0 1 0];
b=[1/6 1/3 1/3 1/6]';
c=[0 1/2 1/2 1]';
[t,y]=RK(fun, tspan, y0, n, A, b, c);
figure (1)
plot(t, y(1,:), 'k-'), hold
plot(t, y(2,:), 'k-.')
legend('resources','population')
figure (2)
plot(y(1,:),y(2,:))
title('orbits')
and this is the 3rd one, called driver_pred_prey. The second one has been changed so now the syntax should be right but yes, when I try to run the third(which calls the first two) basically I get reported this:
• Index exceeds the number of array elements (1).
• Error in RK (line 8)
• k(:,1)=feval(fun,t(j),y(:,j));
• Error in driver_pred_prey (line 14)
• [t,y]=RK(fun, tspan, y0, n, A, b, c);
Those are all the error messages MatLab reports, but I really don't see where is the problem!

Jan on 24 Jul 2019
Index exceeds the number of array elements (1).
Error in RK (line 8)
k(:,1)=feval(fun,t(j),y(:,j));
After the line:
t(1)=tspan(1);
t is a scalar.
If you run this code:
for J=1:n
k(:,1)=feval(fun,t(j),y(:,j));
another error is expected: The for loop uses an uppercase J as counter, but inside the loop there is a lower case j. Because it is not defined, j is the imaginary unit and cannot be used as index.
If you fix this to for j = 1:n, there is the problem, that t is a scalar and you cannot access t(2).
t = linspace(tspan(1), tspan(2), n);

DAVIDE TRONO on 24 Jul 2019
Well as you said there's no problem when defining
t=linspace(tspan(1),tspan(2),n);
but now a new problem is reported:
• Index in position 2 exceeds array bounds (must not exceed 1).
(the other two are the same as before). I guess there's just a little mistake with the indexes.
Jan on 24 Jul 2019
Again: Please post a complete copy of the error message. The details matter.
In which line does the error occur? There is an index in position 2, which exceeds the limits. This should help you to find and to fix the problem.
With some guessing: The code tries to access y(:,j), but this has not been defined anywhere before. The Runge Kutta code calculates the k 's only, but does not create the values for the actual trajectory.
You find many working Runge-Kutta-codes in the net.