there are something wrong with solving DAEs by using ode15s

when i used the ode15s to solve the DAEs,I got some error infomations. here is the equs:

and at the same times x,y,z satisfie an algebraic relation:

so,I got a DAEs system.

here is my codes:

function out =daestest(t,k)

w2=302.*2.*pi./60;
sita=0.79550314./180.*pi;
beta=35.10858727./180.*pi;
w1=(sin(sita+beta)./sin(beta)-1).*w2;
w1_real=w1+w2;
R=995.68198879/1000;
XA=313.33036053/1000;
XB=699.99999999/1000;
ZA=550/1000;
ZB=0;
m=1;
mu=0.01;
g=-9.8;
G=g*m;

a1=(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita));
a2=-(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita));
a3=-cos(t*w2)*sin(sita);
a4=-(-w2^2*(k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))-k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))+R*cos(t*w2)*sin(sita)-k(3)*cos(t*w2)*sin(sita))-2*w2*(k(5)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))-k(4)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))+w1*k(1)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))+w1*k(2)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-k(6)*sin(t*w2)*sin(sita))-2*w1*k(4)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))-2*w1*k(5)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))-w1^2*k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+w1^2*k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))+1/m*((mu*k(7)*k(4)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(5)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(2*k(7)*k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(2*k(7)*k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(mu*k(7)*k(6)*cos(t*w2)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(k(7)*cos(t*w2)*sin(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
b1=-(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita));
b2=(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita));
b3=-sin(t*w2)*sin(sita);
b4=-(-2*w2*(k(5)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))-k(4)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+w1*k(1)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))+w1*k(2)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+k(6)*cos(t*w2)*sin(sita))+w2^2*(k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))-R*sin(t*w2)*sin(sita)+k(3)*sin(t*w2)*sin(sita))+2*w1*k(4)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))+2*w1*k(5)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))+w1^2*k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-w1^2*k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))+1/m*((mu*k(7)*k(5)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(4)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(2*k(7)*k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)+(2*k(7)*k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(mu*k(7)*k(6)*sin(t*w2)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(k(7)*sin(t*w2)*sin(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
c1=cos(t*w1_real)*sin(sita);
c2=-sin(t*w1_real)*sin(sita);
c3=cos(sita);
c4=-(-w1^2*k(1)*cos(t*w1_real)*sin(sita)+w1^2*k(2)*sin(t*w1_real)*sin(sita)-2*w1*k(5)*cos(t*w1_real)*sin(sita)-2*w1*k(4)*sin(t*w1_real)*sin(sita))+1/m*(G+(mu*k(7)*k(6)*cos(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(5)*sin(t*w1_real)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(2*k(7)*k(1)*cos(t*w1_real)*sin(sita))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(2*k(7)*k(2)*sin(t*w1_real)*sin(sita))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)+(mu*k(7)*k(4)*cos(t*w1_real)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(k(7)*cos(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
out=[k(4)
         k(5)
         k(6)
         (a2*b3*c4-a2*b4*c3-a3*b2*c4+a3*b4*c2+a4*b2*c3-a4*b3*c2)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
         (a1*b4*c3-a1*b3*c4+a3*b1*c4-a3*b4*c1-a4*b1*c3+a4*b3*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
         (a1*b2*c4-a1*b4*c2-a2*b1*c4+a2*b4*c1+a4*b1*c2-a4*b2*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
         k(1).^2+k(2).^2-((k(3)-ZB).*(XA-XB)./(ZA-ZB)+XB).^2];

when i used this code to solve the DAEs

y0=[0.3133283972712190;0.001104257865152230;0.5500000247589403;0;0;0;0;];
tspan=[2.533037911852984:5e-4:2.6];
M=[1 0 0 0 0 0 0;0 1 0 0 0 0 0;0 0 1 0 0 0 0;0 0 0 1 0 0 0;0 0 0 0 1 0 0;0 0 0 0 0 1 0;0 0 0 0 0 0 0;];
options = odeset('Mass',M);
[t,k]=ode15s(@daestest,tspan,y0,options);

I got the error:

Error using daeic12 (line 166)
Need a better guess y0 for consistent initial conditions.

Error in ode15s (line 310)
    [y,yp,f0,dfdy,nFE,nPD,Jfac] =
    daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in DAES15s (line 5)
[t,k]=ode15s(@daestest,tspan,y0,options);

So,I changed the initial value of y0 to

y0=[0.3133283972712190;0.001104257865152230;0.5500000247589403;1;1;1;1000;];

I also got the error:

Error using daeic12 (line 76)
This DAE appears to be of index greater than 1.

Error in ode15s (line 310)
    [y,yp,f0,dfdy,nFE,nPD,Jfac] =
    daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in DAES15s (line 5)
[t,k]=ode15s(@daestest,tspan,y0,options);

I'm very confused，What is the cause of such error information? Please help me analyze this problem，thak you.

[Torsten]

7 equations for 6 unknowns ? This cannot work.

Best wishes

Torsten.

[Zilong Zhang]

Hi,thanks for your help.

There are 7 variables：x,y,z,x',y',z' and Fn.

They're all functions of t:x(t),y(t),z(t),x'(t),y'(t),z'(t),Fn(t).

And, in the DAEs functions,the 7 variables are replaced by k(1)to k(7).

Is the error caused by the variable replacement method?

Could you help me analyze it?

Thanks again.

[Torsten]

Try "reduceDAEIndex":

https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html#bvh12tx-2

Best wishes

Torsten.

[Zilong Zhang]

The method in “”is great! I modified my program according this method，here is my code：

syms x(t) y(t) z(t)  m w1 w1_real w2 m R sita Fn(t) G mu XA XB ZA ZB
eqn1=-m.*(2.*w2.*(diff(y(t),1).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))-diff(x(t),1).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))+w1.*x(t).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))+w1.*y(t).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))-diff(z(t),1).*sin(t.*w2).*sin(sita))+w2.^2.*(x(t).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))-y(t).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita))+R.*cos(t.*w2).*sin(sita)-z(t).*cos(t.*w2).*sin(sita))-diff(x(t),2).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))+diff(y(t),2).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita))+2.*w1.*diff(x(t),1).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita))+2.*w1.*diff(y(t),1).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))+diff(z(t),2).*cos(t.*w2).*sin(sita)+w1.^2.*x(t).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))-w1.^2.*y(t).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita)))==(mu.*Fn(t).*diff(x(t),1).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita)))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)-(mu.*Fn(t).*diff(y(t),1).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita)))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)+(2.*Fn(t).*x(t).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita)))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)-(2.*Fn(t).*y(t).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita)))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)-(mu.*Fn(t).*diff(z(t),1).*cos(t.*w2).*sin(sita))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)+(Fn(t).*cos(t.*w2).*sin(sita).*(2.*XB-(2.*(XA-XB).*(ZB-z(t)))./(ZA-ZB)).*(XA-XB))./((ZA-ZB).*(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2));
eqn2=m.*(diff(y(t),2).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))-diff(x(t),2).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))-2.*w2.*(diff(y(t),1).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita))-diff(x(t),1).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))+w1.*x(t).*(cos(t.*w1_real).*sin(t.*w2)-cos(t.*w2).*sin(t.*w1_real).*cos(sita))+w1.*y(t).*(sin(t.*w2).*sin(t.*w1_real)+cos(t.*w2).*cos(t.*w1_real).*cos(sita))+diff(z(t),1).*cos(t.*w2).*sin(sita))+w2.^2.*(x(t).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))-y(t).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))-R.*sin(t.*w2).*sin(sita)+z(t).*sin(t.*w2).*sin(sita))+2.*w1.*diff(x(t),1).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))+2.*w1.*diff(y(t),1).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))+w1.^2.*x(t).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita))-w1.^2.*y(t).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita))-diff(z(t),2).*sin(t.*w2).*sin(sita))==(mu.*Fn(t).*diff(y(t),1).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita)))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)-(mu.*Fn(t).*diff(x(t),1).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita)))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)-(2.*Fn(t).*x(t).*(cos(t.*w2).*sin(t.*w1_real)-cos(t.*w1_real).*sin(t.*w2).*cos(sita)))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)+(2.*Fn(t).*y(t).*(cos(t.*w2).*cos(t.*w1_real)+sin(t.*w2).*sin(t.*w1_real).*cos(sita)))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)-(mu.*Fn(t).*diff(z(t),1).*sin(t.*w2).*sin(sita))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)+(Fn(t).*sin(t.*w2).*sin(sita).*(2.*XB-(2.*(XA-XB).*(ZB-z(t)))./(ZA-ZB)).*(XA-XB))./((ZA-ZB).*(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2));
eqn3=-m.*(diff(y(t),2).*sin(t.*w1_real).*sin(sita)-diff(x(t),2).*cos(t.*w1_real).*sin(sita)-diff(z(t),2).*cos(sita)+w1.^2.*x(t).*cos(t.*w1_real).*sin(sita)-w1.^2.*y(t).*sin(t.*w1_real).*sin(sita)+2.*w1.*diff(y(t),1).*cos(t.*w1_real).*sin(sita)+2.*w1.*diff(x(t),1).*sin(t.*w1_real).*sin(sita))==G+(mu.*Fn(t).*diff(z(t),1).*cos(sita))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)-(mu.*Fn(t).*diff(y(t),1).*sin(t.*w1_real).*sin(sita))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)+(2.*Fn(t).*x(t).*cos(t.*w1_real).*sin(sita))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)-(2.*Fn(t).*y(t).*sin(t.*w1_real).*sin(sita))./(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2)+(mu.*Fn(t).*diff(x(t),1).*cos(t.*w1_real).*sin(sita))./(diff(x(t),1).^2+diff(y(t),1).^2+diff(z(t),1).^2).^(1./2)-(Fn(t).*cos(sita).*(2.*XB-(2.*(XA-XB).*(ZB-z(t)))./(ZA-ZB)).*(XA-XB))./((ZA-ZB).*(4.*x(t).^2+4.*y(t).^2+(4.*(XA-XB).^2.*(XB-((XA-XB).*(ZB-z(t)))./(ZA-ZB)).^2)./(ZA-ZB).^2).^(1./2));
eqn4=x(t)^2+y(t)^2-((z(t)-ZB)*(XA-XB)/(ZA-ZB)+XB)^2==0;
eqns = [eqn1 eqn2 eqn3 eqn4];
vars = [x(t);y(t);z(t);Fn(t)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars)
[eqns, vars] = reduceDifferentialOrder(eqns, vars)
isLowIndexDAE(eqns,vars)
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)
isLowIndexDAE(DAEs,DAEvars)
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars)
f = daeFunction(DAEs, DAEvars, G, R, XA, XB, ZA, ZB, m, mu, sita, w1, w2, w1_real);
m=1;
G=-9.81.*m;
R=995.68198879/1000;
XA=313.33036053/1000;
XB=699.99999999/1000;
ZA=550/1000;
ZB=0;
w2=302.*2.*pi./60;
sita=0.79550314./180.*pi;
beta=35.10858727./180.*pi;
w1=(sin(sita+beta)./sin(beta)-1).*w2;
w1_real=w1+w2;
mu=0.1;
F = @(t, Y, YP) f(t, Y, YP,G, R, XA, XB, ZA, ZB, m, mu, sita, w1, w2, w1_real);
DAEvars
y0est = [0.3133283972712190; 0.001104257865152230; 0.55; 0; 0; 0; 0; 0; 0; 0;];
yp0est = zeros(10,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(F, 2.533037911852984, y0est, [], yp0est, [], opt)

But I got a new error：

Error using decic>sls (line 170)
Index may be greater than one.

Error in decic (line 77)
    [dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);

Error in DAEsL (line 37)
[y0, yp0] = decic(F, 2.533037911852984, y0est, [],
yp0est, [], opt)

why this happend? and what is the meaning of Index？

Torsten,thans for your great help.

[Torsten]

The problem with your system is that there is no explicit ODE equation for k(7), but it must be adjusted such that the algebraic constraint x^2+y^2+...= 0 is satisfied.

My suggestion would be to differentiate the algebraic constraint twice with repsect to t, insert the expressions for x'', y'' and z'' (thus f1(...),f2(...),f3(...)) and solve for k(7). This way, you can circumvent k(7) as solution variable, and it is k(7) which causes the index problem.

Best wishes

Torsten.

[Zilong Zhang]

Thank you Torsten.

In my equations, Fn is a variable about t:Fn=Fn(t).So what can I used to solve Fn?

Because my English is not good, so I can't quite understand the description in "My suggestion would be to differentiate the algebraic constraint twice with repsect to t, insert the expressions for x'', y'' and z'' (thus f1(...),f2(...),f3(...)) and solve for k(7). " can you explain it in more detail?

[Torsten]

Is Fn = k(7) ?

Best wishes

Torsten.

[Zilong Zhang]

Yes,the fn is replaced by k(7). And,x,y,z,x',y' and z' are replaced by k(1),k(2) ,k(3),k(4),k(5) and k(6).

Answer 1

[Torsten]

Differentiating your constraint you get

Differentiating once:

2*x*x'+2*y*y'-2*((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z' = 0

Differentiating twice:

2*x*x''+2*x'^2+2*y*y''+2*y'^2-2*(((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z''+((XA-XB)/(ZA-ZB)*z')^2)=0

Now

x''=f1(x,y,z,x',y',z',Fn,t),
y''=f2(x,y,z,x',y',z',Fn,t) and 
z''=f3(x,y,z,x',y',z',Fn,t) .

Insert these expressions for x'',y'' and z'' in the twice-differentiated constraint condition.

This gives you an algebraic equation for Fn.

Use this algebraic equation to determine Fn instead of

x^2+y^2-((z-ZB)*(XA-XB)/(ZA-ZB)+XB)^2 = 0

in your actual formulation of the DAE-system.

This way, you can be sure that you only work with an index-one problem.

And be careful with the initial conditions:

They have to satisfy the constraint condition as well as the condition differentiated once and twice. Thus you don't have 7 initial conditions you can prescribe freely, but only 7-3=4.

Best wishes

Torsten.

[Zilong Zhang]

Torsten，thanks for your help.

If I use this method,I can get an equation of Fn(x,y,z,z',y',z').

But in the writing rules of DAEs of ODEs,M*y'=Y. The M is mass matrix，the y' are the first differentiate of variables， and the Y are equations of each variables. so,how can I use the equation of Fn(x,y,z,z',y',z')? Should I differentiate the Fn(x,y,z,z',y',z') once, and put the differentiate equations in the right of the "M*y'=Y", just like Fn'=Fn'(x,y,z,z',y',z')?

Thanks very much for your sustained help.

Zilong.

[Torsten]

No, first make an attempt without differentiating the equation for Fn with respect to t.

Just leave the mass matrix you used above as it is and include the equation for Fn as an algebraic equation:

2*x*f1(x,y,z,x',y',z',Fn,t)+2*x'^2+2*y*f2(x,y,z,x',y',z',Fn,t)+2*y'^2-2*...=0

Best wishes

Torsten.

[Zilong Zhang]

Is it correct that use the

2*x*f1(x,y,z,x',y',z',Fn,t)+2*x'^2+2*y*f2(x,y,z,x',y',z',Fn,t)+2*y'^2-2*...=0

to replace the

x^2+y^2-((z-ZB)*(XA-XB)/(ZA-ZB)+XB)^2 = 0

as the forth equations of the DAEs.

And the equations of DAEs system includes:

the mass matrix M of "M*y'=Y":

[1 0 0 0 0 0 0
 0 1 0 0 0 0 0
 0 0 1 0 0 0 0
 0 0 0 1 0 0 0
 0 0 0 0 1 0 0
 0 0 0 0 0 1 0
 0 0 0 0 0 0 0]

the y' of "M*y'=Y":

[x'

   y'

   z'

   x''

   y''

   z''

   Fn]

and the Y of "M*y'=Y":

[k(4)
 k(5)
 k(6)
 eq1
 eq2
 eq3
 2*x*f1(x,y,z,x',y',z',Fn,t)+2*x'^2+2*y*f2(x,y,z,x',y',z',Fn,t)+2*y'^2-2*...]

Is it correct？

[Torsten]

It's not the forth, but the seventh equation of the original DAE system that should be replaced by

2*x*f1(x,y,z,x',y',z',Fn,t)+2*x'^2+2*y*f2(x,y,z,x',y',z',Fn,t)+2*y'^2-2*...

but everything else is correct.

Best wishes

Torsten.

[Zilong Zhang]

Thanks Torsten, it works. But the value of the DAEs are deeply dependent on the initial value of DAEs. And when I use the

   y0=[0.3133283972712190;0.001104257865152230;0.5500000247589403;0.001;0.001;0.001;0;]

as the initial value of DAEs,there are some error:

Error using daeic12 (line 76) This DAE appears to be of index greater than 1.

Error in ode15s (line 310)
    [y,yp,f0,dfdy,nFE,nPD,Jfac] =
    daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in DAES15s (line 6)
[t,k]=ode15s(@daestest,tspan,y0,options);

The reason is caused by using the x,y,z,x',y',z',Fn as the denominators in the forth equation to seventh equation in DAEs? And how do I choose initial values to make the calculation of DAEs more accurate?

here is the latest code:

 the main code:

% y0=[0.3133283972712190;0.001104257865152230;0.5500000247589403;0;0;0;0;];
  y0=[0.3133283972712190;0.001104257865152230;0.5500000247589403;0.001;0.001;0.001;0;];
  tspan=[2.533037911852984:5e-4:2.6];
  M=[1 0 0 0 0 0 0;0 1 0 0 0 0 0;0 0 1 0 0 0 0;0 0 0 1 0 0 0;0 0 0 0 1 0 0;0 0 0 0 0 1 0;0 0 0 0 0 0 0;];
  options = odeset('Mass',M);
  [t,k]=ode15s(@daestest,tspan,y0,options);

    plot3(k(:,1),k(:,2),k(:,3),'b-o')
    grid on
    hold on

the function code:

function out =daestest(t,k)

    w2=302.*2.*pi./60;
    sita=0.79550314./180.*pi;
    beta=35.10858727./180.*pi;
    w1=(sin(sita+beta)./sin(beta)-1).*w2;
    w1_real=w1+w2;
    R=995.68198879/1000;
    XA=313.33036053/1000;
    XB=699.99999999/1000;
    ZA=550/1000;
    ZB=0;
    m=1;
    mu=0.01;
    g=-9.8;
    G=g*m;

    a1=(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita));
    a2=-(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita));
    a3=-cos(t*w2)*sin(sita);
    a4=-(-w2^2*(k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))-k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))+R*cos(t*w2)*sin(sita)-k(3)*cos(t*w2)*sin(sita))-2*w2*(k(5)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))-k(4)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))+w1*k(1)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))+w1*k(2)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-k(6)*sin(t*w2)*sin(sita))-2*w1*k(4)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))-2*w1*k(5)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))-w1^2*k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+w1^2*k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))+1/m*((mu*k(7)*k(4)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(5)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(2*k(7)*k(1)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(2*k(7)*k(2)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(mu*k(7)*k(6)*cos(t*w2)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(k(7)*cos(t*w2)*sin(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
    b1=-(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita));
    b2=(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita));
    b3=-sin(t*w2)*sin(sita);
    b4=-(-2*w2*(k(5)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))-k(4)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+w1*k(1)*(cos(t*w1_real)*sin(t*w2)-cos(t*w2)*sin(t*w1_real)*cos(sita))+w1*k(2)*(sin(t*w2)*sin(t*w1_real)+cos(t*w2)*cos(t*w1_real)*cos(sita))+k(6)*cos(t*w2)*sin(sita))+w2^2*(k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))-R*sin(t*w2)*sin(sita)+k(3)*sin(t*w2)*sin(sita))+2*w1*k(4)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita))+2*w1*k(5)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))+w1^2*k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita))-w1^2*k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))+1/m*((mu*k(7)*k(5)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(4)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita)))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(2*k(7)*k(1)*(cos(t*w2)*sin(t*w1_real)-cos(t*w1_real)*sin(t*w2)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)+(2*k(7)*k(2)*(cos(t*w2)*cos(t*w1_real)+sin(t*w2)*sin(t*w1_real)*cos(sita)))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(mu*k(7)*k(6)*sin(t*w2)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(k(7)*sin(t*w2)*sin(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
    c1=cos(t*w1_real)*sin(sita);
    c2=-sin(t*w1_real)*sin(sita);
    c3=cos(sita);
    c4=-(-w1^2*k(1)*cos(t*w1_real)*sin(sita)+w1^2*k(2)*sin(t*w1_real)*sin(sita)-2*w1*k(5)*cos(t*w1_real)*sin(sita)-2*w1*k(4)*sin(t*w1_real)*sin(sita))+1/m*(G+(mu*k(7)*k(6)*cos(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(mu*k(7)*k(5)*sin(t*w1_real)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)+(2*k(7)*k(1)*cos(t*w1_real)*sin(sita))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)-(2*k(7)*k(2)*sin(t*w1_real)*sin(sita))/(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)+(mu*k(7)*k(4)*cos(t*w1_real)*sin(sita))/(k(4)^2+k(5)^2+k(6)^2)^(1/2)-(k(7)*cos(sita)*(2*XB-(2*(XA-XB)*(ZB-k(3)))/(ZA-ZB))*(XA-XB))/((ZA-ZB)*(4*k(1)^2+4*k(2)^2+(4*(XA-XB)^2*(XB-((XA-XB)*(ZB-k(3)))/(ZA-ZB))^2)/(ZA-ZB)^2)^(1/2)));
    out=[k(4)
             k(5)
             k(6)
             (a2*b3*c4-a2*b4*c3-a3*b2*c4+a3*b4*c2+a4*b2*c3-a4*b3*c2)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
             (a1*b4*c3-a1*b3*c4+a3*b1*c4-a3*b4*c1-a4*b1*c3+a4*b3*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
             (a1*b2*c4-a1*b4*c2-a2*b1*c4+a2*b4*c1+a4*b1*c2-a4*b2*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1)
             2*k(1)*((a2*b3*c4-a2*b4*c3-a3*b2*c4+a3*b4*c2+a4*b2*c3-a4*b3*c2)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1))+2*k(4)^2+2*k(2)*((a1*b4*c3-a1*b3*c4+a3*b1*c4-a3*b4*c1-a4*b1*c3+a4*b3*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1))+2*k(5)^2-2*(((k(3)-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*((a1*b2*c4-a1*b4*c2-a2*b1*c4+a2*b4*c1+a4*b1*c2-a4*b2*c1)/(a1*b2*c3-a1*b3*c2-a2*b1*c3+a2*b3*c1+a3*b1*c2-a3*b2*c1))+((XA-XB)/(ZA-ZB)*k(6))^2)];
    end

[Torsten]

Your initial values must satisfy

x^2+y^2-((z-ZB)*(XA-XB)/(ZA-ZB)+XB)^2 = 0
2*x*x'+2*y*y'-2*((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z' = 0
2*x*x''+2*x'^2+2*y*y''+2*y'^2-2*(((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z''+((XA-XB)/(ZA-ZB)*z')^2) = 0

with x'' = f1(...),y'' = f2(...),z'' = f3(...).

to be reasonable for your problem.

Did you check this ?

Best wishes

Torsten.

[Zilong Zhang]

OK,I'll Try. Thank you very much for your help. These help means a lot to me.

[Torsten]

One further hint:

You should always check whether a solution you obtained remains on your manifold for the complete integration interval:

x^2+y^2-((z-ZB)*(XA-XB)/(ZA-ZB)+XB)^2 = 0

If this is not the case, something is wrong with your code or you have to strengthen the tolerances of ODE15S (RelTol,AbsTol).

Best wishes

Torsten.

[Zilong Zhang]

If i get the initial values of x,y,z, what i can do to find the initial values of x',y',z'to satisfie the equations:

x^2+y^2-((z-ZB)*(XA-XB)/(ZA-ZB)+XB)^2 = 0
2*x*x'+2*y*y'-2*((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z' = 0
2*x*x''+2*x'^2+2*y*y''+2*y'^2-2*(((z-ZB)*(XA-XB)/(ZA-ZB)+XB)*(XA-XB)/(ZA-ZB)*z''+((XA-XB)/(ZA-ZB)*z')^2) = 0

Thanks

[Torsten]

Fix 4 of the 7 unknowns (e.g. x,y,x' and y') and solve the above nonlinear system in the unknowns z,z' and Fn using MATLAB's "fsolve",e.g. .

Best wishes

Torsten.

[Zilong Zhang]

Torsten,you are great,thanks for your help.