How to solve the "Not enough input arguments" in bvp4c ?

Good evening,

I am trying to model the trajectory of a rocket launch for its first stage, from lift-off (0 seconds) to first burn-out (118 seconds), as a boundary value problem using bvp4c. I am using 8 state variables, being y(1)=x, y(2)=x', y(3)=z, y(4)=z', y(5)=v, y(6)=v', y(7)=γ, y(8)=γ'. They make up a system of 8 first order ODEs, so I am using 8 boundary conditions and 8 initial guessed value. When I run the simulation an error comes up, saying "Not enough imput arguments". I have triple checked the code to ensure all variables used appear in the arguments, but it still doesn't solve it.

Here is the code, I am using 4 files .m in total.

ValuesVariables.m (here is where I stored all the variables used):

t=linspace(0,118,119);

FBSL=839100;
FBVac=1020400;
FCSL=792400;
FCVac=985600;
MatB=[(log(0+1)+1) 1 FBSL; (log(180000+1)+1) 1 FBVac];
BatB=rref(MatB); BatB13=BatB(1,3); BatB23=BatB(2,3);
MatC=[(log(0+1)+1) 1 FCSL; (log(180000+1)+1) 1 FCVac];
BatC=rref(MatC); BatC13=BatC(1,3); BatC23=BatC(2,3);

CD= 0.2;

rho_0E=1.225;
rho_41E=0.002927;
ED1=4.7;
MatD1 = [1 1 rho_0E; exp(-ED1) 1 rho_41E];
BatD1=rref(MatD1); BatD113=BatD1(1,3); BatD123=BatD1(2,3); 

dB=2.68;
dC=2.95;
A1=4*pi*(dB/2)^2+pi*(dC/2)^2;

muE=3.98199;
REarth=6371000;

MfB=39600;
MfB=39600;
MfB=39600;
MfFregat=5350*0;
Mp=6560;
MS=25200;
MsB=3810;
MsC=6550;
Ms3=2410;
MsFregat=1000*0;
tBO1=118;
tBO2=280;
MT1E=-4.*((MfB/tBO1).*t)-((MfC/tBO2).*t)+(4.*MfB+MfC+Mf3+MfFregat+4.*MsB+MsC+Ms3+MsFregat+Mp);

bvpTEST.m (here is the bvp4c command is):

clear;
ValuesVariables;

init=bvpinit(linspace(0,118,50),@init_0_118);
sol=bvp4c(@rhs_0_118,@bc_0_118,init);
t=linspace(0,118,119);
BS=deval(sol,t);
plot(t,BS);

rhs_0_118.m (here I defined the ODEs):

function rhs=rhs_0_118(t,y,BatB13,BatB23,BatC13,BatC23,BatD113,BatD123,ED1,CD,A1,muE,REarth,MfB,tBO1,MfC,tBO2,MT1E)

rhs=zeros(8,1);

rhs(1)=y(5).*cos(y(7));
rhs(2)=((2.*y(6).*((((y(2)).^2)+((y(4)).^2)).^(0.5)))-(2.*y(4).*y(6).*sin(y(7)))).*((2.*y(2)).^(-1));
rhs(3)=y(5).*sin(y(7));
rhs(4)=((2.*y(6).*((((y(2)).^2)+((y(4)).^2)).^(0.5)))-(2.*y(2).*y(6).*cos(y(7)))).*((2.*y(4)).^(-1));
rhs(5)=(((((BatB13.*(log(y(3)+1)+1)+BatB23))+((BatC13.*(log(y(3)+1)+1)+BatC23)))-((0.5).*CD.*A1.*((y(5)).^2).*(BatD113.*exp(-ED1.*(y(3)).*(1/41000))+BatD123))).*((MT1E).^(-1)))-((muE.*((REarth+y(3)).^(-2))).*cos(y(7)));
rhs(6)=-(((BatB13.*(log(y(3)+1)+1)+BatB23))+((BatC13.*(log(y(3)+1)+1)+BatC23))).*((MT1E).^(-2)).*(-4.*(MfB/tBO1)-(MfC/tBO2))-(0.5).*(BatD113.*exp(-ED1.*(y(3)).*(1/41000))+BatD123).*(CD).*(A1).*(2.*(y(5).*(y(6)).*((MT1E).^(-1)))-((y(5)).^(2)).*((MT1E).^(-2)).*(-4.*(MfB/tBO1)-(MfC/tBO2)))+muE.*((REarth+y(3)).^(-2)).*y(8).*sin(y(7));
rhs(7)=-(muE.*((REarth+y(3)).^(-2))).*(cos(y(7)).*((y(5)).^(-1)));
rhs(8)=(muE.*((REarth+y(3)).^(-2))).*(((y(5)).^(-1)).*y(8).*sin(y(7))+((y(5)).^(-2)).*y(6));

end

bc_0_118.m (here I defined all the boundary conditions):

function bc=bc_0_118(yl,yr)

bc=zeros(8,1);

bc(1)=yl(1);
bc(2)=yl(3);
bc(3)=yl(5);
bc(4)=yl(7)-(pi/2);
bc(5)=yr(1)-39000;
bc(6)=yr(3)-41000;
bc(7)=yr(5)-(8300.*1000.*((60*60).^(-1)));
bc(8)=yr(7)-(pi/3);

end

init_0_118.m (here I defined the initial guessed values):

function in=init_0_118(t)

in=zeros(8,1);

in(1)=1;
in(2)=1;
in(3)=1;
in(4)=9.81;
in(5)=1;
in(6)=9.81;
in(7)=(pi/2);
in(8)=0;

end

That is the whole code, and here is the error:

Not enough input arguments.

Error in rhs_0_118 (line 9)
rhs(5)=(((((BatB13.*(log(y(3)+1)+1)+BatB23))+((BatC13.*(log(y(3)+1)+1)+BatC23)))-((0.5).*CD.*A1.*((y(5)).^2).*(BatD113.*exp(-ED1.*(y(3)).*(1/41000))+BatD123))).*((MT1E).^(-1)))-((muE.*((REarth+y(3)).^(-2))).*cos(y(7)));

Error in bvparguments (line 105)
    testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 130)
    bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in bvpTEST (line 5)
sol=bvp4c(@rhs_0_118,@bc_0_118,init);

I tried including all variable names in the all function files, or only writing the ones that appear in the files assuming that each variable is accounted for since I call the "ValuesVariables" anyway, but it doesn't work. I suspect I didn't write the variable "t" where I should have, and that the problem has to do with the "init_0_118.m" file, but I can't spot it.

Any help with this would be highly appreciated, thank you in advance.

Dave

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Dave B]

Thank you for the quick reply. I had a look at the link you provided and I changed the code. Now I am calling the "ValuesVariables" directly from the rhs_0_118.m function file, which now looks like this:

function rhs=rhs_0_118(t,y,BatB13,BatB23,BatC13,BatC23,BatD113,BatD123,ED1,CD,A1,muE,REarth,MfB,tBO1,MfC,tBO2,Mf3,MfFregat,MsB,MsC,Ms3,MsFregat,Mp,FBSL,FBVac,FCSL,FCVac,MatB,BatB,MatC,BatC,rho_0E,rho_41E,MatD1,BatD1,dB,dC)
ValueVariables;

rhs=zeros(8,1);

rhs(1)=y(5).*cos(y(7));
rhs(2)=((2.*y(6).*((((y(2)).^2)+((y(4)).^2)).^(0.5)))-(2.*y(4).*y(6).*sin(y(7)))).*((2.*y(2)).^(-1));
rhs(3)=y(5).*sin(y(7));
rhs(4)=((2.*y(6).*((((y(2)).^2)+((y(4)).^2)).^(0.5)))-(2.*y(2).*y(6).*cos(y(7)))).*((2.*y(4)).^(-1));
       
    velocitydot=@(t)(((BatB13.*(log(y(3)+1)+1)+BatB23)+(BatC13.*(log(y(3)+1)+1)+BatC23))-((0.5).*CD.*A1.*((y(5)).^2).*(BatD113.*exp(-ED1.*(y(3)).*(1/41000))+BatD123))).*((-4.*((MfB/tBO1).*t)-((MfC/tBO2).*t)+(4.*MfB+MfC+Mf3+MfFregat+4.*MsB+MsC+Ms3+MsFregat+Mp)).^(-1))-(muE.*((REarth+y(3)).^(-2)).*cos(y(7)));

rhs(5)=velocitydot(t);

    velocity2dots=@(t)-(((BatB13.*(log(y(3)+1)+1)+BatB23))+((BatC13.*(log(y(3)+1)+1)+BatC23))).*((-4.*((MfB/tBO1).*t)-((MfC/tBO2).*t)+(4.*MfB+MfC+Mf3+MfFregat+4.*MsB+MsC+Ms3+MsFregat+Mp)).^(-2)).*(-4.*(MfB/tBO1)-(MfC/tBO2))-(0.5).*(BatD113.*exp(-ED1.*(y(3)).*(1/41000))+BatD123).*(CD).*(A1).*(2.*(y(5).*(y(6)).*((-4.*((MfB/tBO1).*t)-((MfC/tBO2).*t)+(4.*MfB+MfC+Mf3+MfFregat+4.*MsB+MsC+Ms3+MsFregat+Mp)).^(-1)))-((y(5)).^(2)).*((-4.*((MfB/tBO1).*t)-((MfC/tBO2).*t)+(4.*MfB+MfC+Mf3+MfFregat+4.*MsB+MsC+Ms3+MsFregat+Mp)).^(-2)).*(-4.*(MfB/tBO1)-(MfC/tBO2)))+muE.*((REarth+y(3)).^(-2)).*y(8).*sin(y(7));
    
rhs(6)=velocity2dots(t);
rhs(7)=-(muE.*((REarth+y(3)).^(-2))).*(cos(y(7)).*((y(5)).^(-1)));
rhs(8)=(muE.*((REarth+y(3)).^(-2))).*(((y(5)).^(-1)).*y(8).*sin(y(7))+((y(5)).^(-2)).*y(6));

end

However, it seems that it still does not progress due to the following:

Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in bvpTEST (line 5)
sol=bvp4c(@rhs_0_118,@bc_0_118,init);

Is there a chance I interpreted the link wrong?

-Dave

[Torsten]

sol=bvp4c(@(t,y)rhs_0_118(t,y,BatB13,BatB23,BatC13,BatC23,BatD113,BatD123,...
                ED1,CD,A1,muE,REarth,MfB,tBO1,MfC,tBO2,Mf3,MfFregat,MsB,MsC,...
                Ms3,MsFregat,Mp,FBSL,FBVac,FCSL,FCVac,MatB,BatB,MatC,BatC,...
                rho_0E,rho_41E,MatD1,BatD1,dB,dC),@bc_0_118,init);
            

And why do you create function handles "velocitydot" and "velocity2dots" ? Just insert the expressions for rhs(5) and rhs(6).

[Dave B]

Hi, thank you for the tip. I tried writing it how you suggested, and even writing the rhs(5) and rhs(6) expressions directly but the single Jacobian problem persists. Could there be a problem in the way Matlab handles my boundary and initial conditions?

[Torsten]

Do you get reasonable values for rhs(1) - rhs(8) for your initial conditions ? Could you check in rhs_0_118 ?

[Dave B]

I am not sure what the actual value of rhs(8) should be, but plugging in the initial conditions yields a finite (very small) value. rhs(1) gives me zero, and it is what I would expect for the first duration of the rocket launch. rhs(1) is the horizontal velocity, and it should be zero considering that the first stage starts vertically.

I also dubled checked that my boundary and initial conditions don't lead to any denominator being zero, so I changed bc(3)=yl(5); to bc(3)=yl(5)-1; since y(5) appears in some denominators. I also tried not having any zero values at all for boundary or initial conditions, but I still get the singular Jacobian.