there are something wrong with solving DAEs by using ode15s

2017 11 23
1 个回答
回答已采纳
更新时间：2017 11 28
125 次查看（30 天）

1 个投票

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.

8 个评论
显示 6更早的评论 隐藏 6更早的评论

Torsten
Torsten 2017-11-23
7 equations for 6 unknowns ? This cannot work.
Best wishes
Torsten.
Zilong Zhang
Zilong Zhang 2017-11-23
编辑：Zilong Zhang 2017-11-23
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
Torsten 2017-11-23
Try "reduceDAEIndex":
https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html#bvh12tx-2
Best wishes
Torsten.
Zilong Zhang
Zilong Zhang 2017-11-23
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
Torsten 2017-11-23
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
Zilong Zhang 2017-11-23
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
Torsten 2017-11-23
Is Fn = k(7) ?
Best wishes
Torsten.
Zilong Zhang
Zilong Zhang 2017-11-23
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).

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2017-11-24

0 个投票

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.

11 个评论
显示 9更早的评论 隐藏 9更早的评论

Zilong Zhang
Zilong Zhang 2017-11-24
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
Torsten 2017-11-24
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
Zilong Zhang 2017-11-24
编辑：Zilong Zhang 2017-11-24
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
Torsten 2017-11-24
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
Zilong Zhang 2017-11-24
编辑：Zilong Zhang 2017-11-24
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
Torsten 2017-11-24
编辑：Torsten 2017-11-24
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
Zilong Zhang 2017-11-24
OK,I'll Try. Thank you very much for your help. These help means a lot to me.
Torsten
Torsten 2017-11-24
编辑：Torsten 2017-11-24
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
Zilong Zhang 2017-11-26
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
Torsten 2017-11-27
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
Zilong Zhang 2017-11-28
Torsten,you are great,thanks for your help.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 配置仿真条件 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!