Implementation of RK45 (system of three equations(​dx/dt,dy/d​t,dz/dt))(​errors in imput arguments and number of array elements)

Greetings, I am working on a Runge-Kutta-Fehlberg Method (RKF45) for a system of three equations with respect to time and really need some help with some errors (dx = a * and * xa * (c * x + (1./2) * (dc) * ( abs (x + 1) -abs (x-1))); dy = xyz; dz = -b * y;), the error is "Not enough input arguments. Error in Rk45B> @ (x, y, z) xyz (line 53) fy = @ (x, y, z) xyz; Error in Rk45B (line 116) k1y = h * feval (fy, T (j), Y (j)); ", to verify I change the equation dy, and the error becomes "Index exceeds the number of array elements (1). Error in Rk45B (line 112) k1x = h * feval (fx, T (j), X (j));" , I also make sure to try to correct that error by replacing j by j + 1 and the error persists. the method is based on the RK45 of the book "Numerical Methods Using MATLAB Fourth Edition, page 511". here is my code with the j+1 fix that i made but with the original dy equation.

function[X , Y , Z, T]=Rk45B(xa,ya,za,M,tol) 
%caracteristic of the program
%Input - f is the function entered as a string ’f’ 
% - a and b are the left and right endpoints 
% - ya is the initial condition y(a) 
% - M is the number of steps 
% - tol is the tolerance 
%Output - R=[T’ Y’] where T is the vector of abscissas 
% and Y is the vector of ordinates 
%Enter the coefficients necessary to calculate the 
%values in (28) and (29)
%Rk45B(-1.6,0,1.6,15,25.58,-0.7142857,-1.1428571,100,0.0000000001);



%% this way doesnt work
%function fy = myfx(x,y,z)
%l=-1.6;
%r=0;
%c=1.6;
%d=15;
%  f = @(x,y,z) l*y-l*x-l*(c*x+(1./2)*(d-c)*(abs(x+1)-abs(x-1)));
%  fx = f(x,y,z);
%end
%function fx = myfy(x,y,z)
    %l=-1.6;
    %r=0;
    %c=1.6;
    %d=15; 
 %  f = @(x,y,z) x-y-z;
 %  fy = f(x,y,z);
%end
%function fz = myfz(x,y,z)
    %l=-1.6;
    %r=0;
    %c=1.6;
    %d=15;
 %  f = @(x,y,z) -r*y;
 %  fz = f(x,y,z);
%end
%setup_callback(@myfx);
%setup_callback(@myfy);
%setup_callback(@myfz);
%setup_callback(@fx);
%setup_callback(@fy);
%setup_callback(@fz);
%evaluate_callback();
%%
%%In Main run the line [x,y,z,c]=Rk45B(-1.6 ,0 ,1.6,100 ,0.0000000001 ); 
%%

fx = @(x,y,z) 15*y-15*x-15*((-5./7)*x+(1./2)*((-8./7)-(-5./7))*(abs(x+1)-abs(x-1)));
fy = @(x,y,z) x-y-z;
fz = @(x,y,z) -25.58*y;  

a=0;
b=100;

a2=1/4;
b2=1/4;
a3=3/8;
b3=3/32;
c3=9/32;
a4=12/13; 
b4=1932/2197;
c4=-7200/2197;
d4=7296/2197;
a5=1;
b5=439/216;
c5=-8;
d5=3680/513;
e5=-845/4104;
a6=1/2;
b6=-8/27;
c6=2;
d6=-3544/2565;
e6=1859/4104;
f6=-11/40;
r1=1/360;
r3=-128/4275;
r4=-2197/75240;
r5=1/50;
r6=2/55;
n1=25/216;
n3=1408/2565;
n4=2197/4104;
n5=-1/5;

big=1e15; 
h=(b-a)/M;
hmin=h/64;
hmax=64*h; 
max1=200;

X(1)=xa;
Y(1)=ya;
Z(1)=za;

T(1)=a;
j=1;
br=b-0.00001*abs(b);


while (T(j)<b) 
    
    if ((T(j)+h)>br) 
        h=b-T(j); 
    end
    
    
    %% previous state   T(j),X(j) same for the other kn equations  
    k1x=h*feval(fx,T(j+1),X(j+1)); 
    x2=X(j)+b2*k1x;
    if big<abs(x2)break,end
   %% previous state   T(j),Y(j) 
    k1y=h*feval(fy,T(j+1),Y(j+1)); 
    y2=Y(j)+b2*k1y;
    if big<abs(y2)break,end
    %% previous state   T(j),Z(j) 
    k1z=h*feval(fz,T(j+1),Z(j+1)); 
    z2=Z(j)+b2*k1z;
    if big<abs(z2)break,end    
    
    
    %%
    k2y=h*feval(fy,T(j+1)+a2*h,y2);
    y3=Y(j)+b3*k1y+c3*k2y;
    if big<abs(y3)break,end

    k2x=h*feval(fx,T(j+1)+a2*h,x2);
    x3=X(j)+b3*k1x+c3*k2x;
    if big<abs(x3)break,end
    
    k2z=h*feval(fz,T(j+1)+a2*h,z2);
    z3=Z(j)+b3*k1z+c3*k2z;
    if big<abs(z3)break,end    
    
    %%
    k3y=h*feval(fy,T(j+1)+a3*h,y3);
    y4=Y(j)+b4*k1y+c4*k2y+d4*k3y;
    if big<abs(y4)break,end

    k3x=h*feval(fx,T(j+1)+a3*h,x3);
    x4=X(j)+b4*k1x+c4*k2x+d4*k3x;
    if big<abs(x4)break,end    

    k3z=h*feval(fz,T(j+1)+a3*h,z3);
    z4=Z(j)+b4*k1z+c4*k2z+d4*k3z;
    if big<abs(z4)break,end    
    
    %%
    k4y=h*feval(fy,Y(j+1)+a4*h,y4);
    y5=Y(j)+b5*k1y+c5*k2y+d5*k3y+e5*k4y;
    if big<abs(y5)break,end
    
    k4x=h*feval(fx,X(j+1)+a4*h,x4);
    x5=X(j)+b5*k1x+c5*k2x+d5*k3x+e5*k4x;
    if big<abs(x5)break,end
    
    k4z=h*feval(fz,Z(j+1)+a4*h,z4);
    z5=Z(j)+b5*k1z+c5*k2z+d5*k3z+e5*k4z;
    if big<abs(z5)break,end    
   
    %%
    k5y=h*feval(fy,T(j+1)+a5*h,y5);
    y6=Y(j)+b6*k1y+c6*k2y+d6*k3y+e6*k4y+f6*k5y;
    if big<abs(y6)break,end
    
    k5x=h*feval(fx,T(j+1)+a5*h,x5);
    x6=X(j)+b6*k1x+c6*k2x+d6*k3x+e6*k4x+f6*k5x;
    if big<abs(x6)break,end
    
    k5z=h*feval(fz,T(j+1)+a5*h,z5);
    z6=Z(j)+b6*k1z+c6*k2z+d6*k3z+e6*k4z+f6*k5z;
    if big<abs(z6)break,end    
    
    
    
    %%
    k6y=h*feval(fy,Y(j)+a6*h,y6);
    k6x=h*feval(fx,X(j)+a6*h,x6);
    k6z=h*feval(fz,Z(j)+a6*h,z6);
    
    %%
    erry=abs(r1*k1y+r3*k3y+r4*k4y+r5*k5y+r6*k6y);
    errx=abs(r1*k1x+r3*k3x+r4*k4x+r5*k5x+r6*k6x);
    errz=abs(r1*k1z+r3*k3z+r4*k4z+r5*k5z+r6*k6z);
    
    
    ynew=Y(j)+n1*k1y+n3*k3y+n4*k4y+n5*k5y; 
    xnew=X(j)+n1*k1x+n3*k3x+n4*k4x+n5*k5x;
    znew=Z(j)+n1*k1z+n3*k3z+n4*k4z+n5*k5z;
    
    
   %%
   if (erry<=errx && erry<=errz)
       err =erry;
       A=Y;
       new=ynew;
     else
        if(errx<erry && errx<errz)
           err =errx;
           A=X;
           new=xnew;
          else
           err =errz;
           A=Z;
           new=znew;
        end    
   end
   
   
   %% A , err y new debe ser la variable escogida en el intervalo digase X Y o Z
   
   
   
    %%
    %Error and step size control 
    if((err<tol)||(h<2*hmin)) 
        
      A(j+1)=new; 
      
      if((T(j)+h)>br) 
    
         T(j+1)=b; 
    
      else
         T(j+1)=T(j)+h; 
      end
      j=j+1;
    end
    
if (err==0)     
 s=0;    
 else
    s=0.84*(tol*h/err)^(0.25);
end

if((s<0.75)&&(h>2*hmin)) 
    
   h=h/2; 
end

if((s>1.50)&&(2*h<hmax)) 
   h=2*h; 
end


if((big<abs(A(j)))||(max1==j)),break,end 
M=j; 
if (b>T(j)) 
   M=j+1; 
else
    M=j+1;
end

end
%%

 disp(X(j));
 disp(Y(j));
 disp(Z(j));


end

[Walter Roberson]

You define fx, fy, fz as each needing 3 parameters. When you feval() them, you are supplying only two parameters.

Note: if you have something that is known to be a function handle, instead of using feval(), just invoke the function handle. For example, instead of

feval(fx, X(j)+a6*h, y6)

invoke

fx(X(j)+a6*h, y6, z6)

[darova]

Why the number of input arguments is 2? All you function fx, fy, fz have 3

feval(fz,T(j+1)+a2*h,z2)    % passing time instead of X?

At the end of loop you have:

ynew=Y(j)+n1*k1y+n3*k3y+n4*k4y+n5*k5y; 

But it never used

[Jayser Blanco]

ok ok , i get it. i will fix the inputs arguments in the code ,with the right invoke. ynew is used in the if that follow that line, its used to calculate the smaller step zise in each iteration but i will check that to be more careful.I leave that definition(time instead of Z) because time is part of the theoretical formulas of the RK45, but i think that its no neccesary in that part.

[Jayser Blanco]

ok i fix and all good , but now im getting "Conversion to cell from double is not possible. Error in Rk45B (line 118) X(1)=xa;", but if i erase the prellocations that matlab suggest the program run but take a lot of time.

function[X , Y , Z, T]=Rk45B(xa,ya,za,M,tol) 
%caracteristic of the program
%Input - f is the function entered as a string ’f’ 
% - a and b are the left and right endpoints 
% - ya is the initial condition y(a) 
% - M is the number of steps 
% - tol is the tolerance 
%Output - R=[T’ Y’] where T is the vector of abscissas 
% and Y is the vector of ordinates 
%Enter the coefficients necessary to calculate the 
%values in (28) and (29)
%Rk45B[x,y,z,c]=Rk45B(-1.6 ,0 ,1.6,10 ,0.000001 );


fx = @(x,y,z) 15*y-15*x-15*((-5./7)*x+(1./2)*((-8./7)-(-5./7))*(abs(x+1)-abs(x-1)));
fy = @(x,y,z) x-y-z;
fz = @(x,y,z) -25.58*y;  

a=0;
b=100;

a2=1/4;
b2=1/4;
a3=3/8;
b3=3/32;
c3=9/32;
a4=12/13; 
b4=1932/2197;
c4=-7200/2197;
d4=7296/2197;
a5=1;
b5=439/216;
c5=-8;
d5=3680/513;
e5=-845/4104;
a6=1/2;
b6=-8/27;
c6=2;
d6=-3544/2565;
e6=1859/4104;
f6=-11/40;
r1=1/360;
r3=-128/4275;
r4=-2197/75240;
r5=1/50;
r6=2/55;
n1=25/216;
n3=1408/2565;
n4=2197/4104;
n5=-1/5;

big=1e15; 
h=(b-a)/M;
hmin=h/64;
hmax=64*h; 
max1=200;

X= cell(1,9000);
%X=[];
Y= cell(1,9000);
%Y=[];
Z= cell(1,9000);
%Z=[];
k1x= cell(1,9000);
%k1x=[];
k1y= cell(1,9000);
%k1y=[];
k1z= cell(1,9000);
%k1z=[];
x2= cell(1,9000);
%x2=[];
y2= cell(1,9000);
%y2=[];
z2= cell(1,9000);
%z2=[];
T= cell(1,9000);
%T=[];



X(2)=xa;
Y(2)=ya;
Z(2)=za;

T(1)=a;
j=1;
br=b-0.00001*abs(b);


while (T(j)<b) 
    
    if ((T(j)+h)>br) 
        h=b-T(j); 
    end
    
    
    %%     
    k1x(j+1)=h*fx(X(j+1),Y(j+1),Z(j+1)); 
    x2(j+1)=X(j+1)+b2*k1x(j+1);
    
    if big<abs(x2)
        break
    end

    k1y(j+1)=h*fy(X(j+1),Y(j+1),Z(j+1)); 
    y2(j+1)=Y(j+1)+b2*k1y(j+1);
    if big<abs(y2)
        break
    end
    
    k1z(j+1)=h*fz(X(j+1),Y(j+1),Z(j+1)); 
    z2(j+1)=Z(j+1)+b2*k1z(j+1);
    if big<abs(z2)
        break
    end    
    
    
    %%
    k2y=h*fy(x2,y2,z2);
    y3=Y(j+1)+b3*k1y+c3*k2y;
    if big<abs(y3)
        break
    end

    k2x=h*fx(x2,y2,z2);
    x3=X(j+1)+b3*k1x+c3*k2x;
    if big<abs(x3)
        break
    end
    
    k2z=h*fz(x2,y2,z2);
    z3=Z(j+1)+b3*k1z+c3*k2z;
    if big<abs(z3)
        break
    end    
    
    %%
    k3y=h*fy(x3,y3,z3);
    y4=Y(j+1)+b4*k1y+c4*k2y+d4*k3y;
    if big<abs(y4)
        break
    end

    k3x=h*fx(x3,y3,z3);
    x4=X(j+1)+b4*k1x+c4*k2x+d4*k3x;
    if big<abs(x4)
        break
    end    

    k3z=h*fz(x3,y3,z3);
    z4=Z(j+1)+b4*k1z+c4*k2z+d4*k3z;
    if big<abs(z4)
        break
    end    
    
    %%
    k4y=h*fy(x4,y4,z4);
    y5=Y(j+1)+b5*k1y+c5*k2y+d5*k3y+e5*k4y;
    if big<abs(y5)
        break
    end
    
    k4x=h*fx(x4,y4,z4);
    x5=X(j+1)+b5*k1x+c5*k2x+d5*k3x+e5*k4x;
    if big<abs(x5)
        break
    end
    
    k4z=h*fz(x4,y4,z4);
    z5=Z(j+1)+b5*k1z+c5*k2z+d5*k3z+e5*k4z;
    if big<abs(z5)
        break
    end    
   
    %%
    k5y=h*fy(x5,y5,z5);
    y6=Y(j+1)+b6*k1y+c6*k2y+d6*k3y+e6*k4y+f6*k5y;
    if big<abs(y6)
        break
    end
    
    k5x=h*fx(x5,y5,z5);
    x6=X(j+1)+b6*k1x+c6*k2x+d6*k3x+e6*k4x+f6*k5x;
    if big<abs(x6)
        break
    end
    
    k5z=h*fz(x5,y5,z5);
    z6=Z(j+1)+b6*k1z+c6*k2z+d6*k3z+e6*k4z+f6*k5z;
    if big<abs(z6)
        break
    end    
    
    
    
    %%
    k6y=h*fy(x6,y6,z6);
    k6x=h*fx(x6,y6,z6);
    k6z=h*fz(x6,y6,z6);
    
    %%
    erry=abs(r1*k1y+r3*k3y+r4*k4y+r5*k5y+r6*k6y);
    errx=abs(r1*k1x+r3*k3x+r4*k4x+r5*k5x+r6*k6x);
    errz=abs(r1*k1z+r3*k3z+r4*k4z+r5*k5z+r6*k6z);
    
    
    ynew=Y(j+1)+n1*k1y+n3*k3y+n4*k4y+n5*k5y; 
    xnew=X(j+1)+n1*k1x+n3*k3x+n4*k4x+n5*k5x;
    znew=Z(j+1)+n1*k1z+n3*k3z+n4*k4z+n5*k5z;
    
    
   %%
   if (erry<=errx & erry<=errz)
       err =erry;
       A=Y;
       new=ynew;
     else
        if(errx<erry & errx<errz)
           err =errx;
           A=X;
           new=xnew;
          else
           err =errz;
           A=Z;
           new=znew;
        end    
   end
   
   
   %% A , err y new debe ser la variable escogida en el intervalo digase X Y o Z
   
   
   
    %%
    %Error and step size control 
    if((err<tol)|(h<2*hmin)) 
        
      A(j+1)=new; 
      
      if((T(j)+h)>br) 
    
         T(j+1)=b; 
    
      else
         T(j+1)=T(j)+h; 
      end
      j=j+1;
    end
    
if (err==0)     
 s=0;    
 else
    s=0.84*(tol*h./err).^(0.25);
end

if((s<0.75)&(h>2*hmin)) 
    
   h=h/2; 
end

if((s>1.50)&(2*h<hmax)) 
   h=2*h; 
end


if((big<abs(A(j)))|(max1==j)),break,end 
M=j; 
if (b>T(j)) 
   M=j+1; 
else
    M=j+1;
end

end
%%

 disp(X(j));
 disp(Y(j));
 disp(Z(j));
 disp(T(j))

end

[Walter Roberson]

X{2} = xa;

or

X(2) = {xa};

but not

X(2) = xa;    %will not work when X is a cell and xa is numeric

[Jayser Blanco]

thanks Walter that fix the convertion to cell problem , the thing now will be how i use the operators "<" and ">" without getting "Undefined operator '<' for input arguments of type 'cell'. while (T(j)<b) " .