I do not understand why my code still gives me errors

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 $$

[infinity]

Hello, there are many errors to fix. If you use bebug in Matlab, you will find out several of them.

For example, in function 2.

function dy=pred_prey(t,y)
             function s=cos(x)

Matlab doesn't understand what is "x" and so on.

[Jan]

I tried to clean up your code. Please use the tools for formatting by your own. A standard indentation is useful also, because it makes the reading easier.

function s=cos(x)
r(t)=int(s(x^2),x,0,t);

Here the "function" is invalid. Maybe you mean a symbolic variable? By the way, the integral of the cos() is easy. Do not let it compute in each iteration.

Whenever you mention, that there is an error message, post a complete copy of it. It is easier to solve a problem than to guess, what theproblem is.

It is unlikely, that the code shows 2 errors, because Matlab must stop after the first one.

[Adam]

We aren't likely to understand your errors either if you tell us nothing about what they are! If they are errors they come with useful error information in the command window.

As an aside, you really shouldn't be having instructions like

close all

in an arbitrary function. From a software design perspective a function whose job it is to do a calculation, even if it is also creating plots has no business just arbitrarily closing all figures, whether they are related to it or not! Especially since in your case you create new figures anyway so it's not as though your function just plots to whatever happens to be the current figure on entering the function.

[DAVIDE TRONO]

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]

What about:

% 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);

Please post the complete messages.

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]

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!

Answer 1

[Jan]

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).

Maybe you want this instead:

      t = linspace(tspan(1), tspan(2), n);

[DAVIDE TRONO]

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]

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.