高手留步，MATLAB解DAE问题

千山之巅

解上图方程，以下是我的求解程序


k=1.1;
n=0.8;
odefun=@(t,x)[-n*pi*x(2)^2*x(4);%t对应theta，x(1)对应V，x（2）对应R，x（3）对应z，想（4）对应dz
    -1/(tan(k*t))*x(4);
    x(4);
    x(1)*(sin(t))^3-pi/3*(x(2)^3)*((2-3*cos(t)+(cos(t))^3))];
tspan=[0.1, pi/2];
x0=[pi/3,1,1,0];
M=eye(4,4);
M(4,4)=0;
options=odeset('Mass',M);
[t,x]=ode23t(odefun,tspan,x0,options);
但总是提示
Error using daeic12 (line 77)This DAE appears to be of index greater than 1.
Error in ode23t (line 289)    [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in bingye (line 12)[t,x]=ode23t(odefun,tspan,x0,options);

TBE_Legend

不会matlab，用mathematica做了下

TBE_Legend

1. t0 = 0.1; tend = Pi/2;
   
2. Manipulate[
   
3. {rsol, zsol, vsol} =
   
4.   NDSolveValue[{r'[t] == (-z'[t])/Tan[ k t],
   
5.      v'[t] == n (Pi r[t]^2 z'[t]),
   
6.      v[t]/r[t]^3 == (Pi/3) ((2 - 3 Cos[t] + Cos[t]^3)/Sin[t]^3),
   
7.      v[t0] == Pi/3, r[t0] == 1, z[t0] == 1}, {r, z, v}, {t, t0,
   
8.      tend}] // Quiet;
   
9. Plot[{rsol[t], zsol[t], vsol[t]}, {t, t0, tend},
   
10.   AxesLabel -> Automatic, PlotTheme -> "Detailed",
    
11.   PlotLegends -> "Expressions"], {k, 0.1, 1}, {n, 0.1, 1}]

复制代码