Problem of Differential Algebraic Equations in Matlab

Question

0 投票

Hi all,

I have written a program in order to solve for 128 differential equations and 64 algebraic equations. I use the ode23t, and set Mass matrix M=[eye(128,192);zeros(64,192)];But error reports are returned as

SWITCH expression must be a scalar or string constant.

Error in odemass (line 66) switch(mass) Please help me with this, thank you very much!! The following is firstly the m-file creating DAEs, then the main function

function dy=ringsystemmanu(t,y)
%set up deltaeq and gain P(i)
alpha=10;
for i=1:64
    deltaeq(i)=2*pi*(i-1)/64;
end
P(1)=2-cos(deltaeq(1)-deltaeq(2))-cos(deltaeq(1)-deltaeq(64))+6*(sin(deltaeq(1)-deltaeq(2))+sin(deltaeq(1)-deltaeq(64)));
P(64)=2-cos(deltaeq(64)-deltaeq(1))-cos(deltaeq(64)-deltaeq(63))+6*(sin(deltaeq(64)-deltaeq(1))+sin(deltaeq(64)-deltaeq(63)));
for j=2:63
    P(j)=2-cos(deltaeq(j)-deltaeq(j+1))-cos(deltaeq(j)-deltaeq(j-1))+6*(sin(deltaeq(j)-deltaeq(j+1))+sin(deltaeq(j)-deltaeq(j-1)));
end
%set up 128 differential equations
h=2;
dy=zeros(192,1);
dy(1)=y(2);
dy(2)=P(1)-6*alpha*sin(y(1)-y(129));
for k=3:2:125
    dy(k)=y(k+1);
    dy(k+1)=P(h)-6*alpha*sin(y(k)-y(h+128));
    h=h+1;
end
dy(127)=y(128);
dy(128)=P(64)-6*alpha*sin(y(127)-y(192));
dy(129)=(2-cos(y(129)-y(130))-cos(y(129)-y(192)))+6*(sin(y(129)-y(130))+sin(y(129)-y(192)))-6*alpha*sin(y(1)-y(129));
for n=130:191
    dy(n)=(2-cos(y(n)-y(n+1))-cos(y(n)-y(n-1)))+6*(sin(y(n)-y(n+1))+sin(y(n)-y(n-1)))-6*alpha*sin(y(2*(n-128)-1)-y(n));
end
dy(192)=2-cos(y(192)-y(129))-cos(y(192)-y(191))+6*(sin(y(192)-y(129))+sin(y(192)-y(191)))-6*alpha*sin(y(127)-y(192));
main function
tspan=[0 2.39 5.22 8];
y0=zeros(192,1);
m=1;
for k=1:2:127
    y0(k)=0.5*exp(-0.1*(m-31.5).^2);
    m=m+1;
end
for n=2:2:128
    y0(n)=0;
end
for r=129:192
    y0(r)=2*pi*(r-128-1)/64;
end
M=[eye(128,192);zeros(64,192)];
options=odeset('mass',M);
[t,y]=ode23t('ringsystemmanu',tspan,y0,options);
data=[t,y];
save ringsystem_data.txt data -ascii
subplot(2,2,1)
plot(y(1,1:2:127))
subplot(2,2,2)
plot(y(2,1:2:127))
subplot(2,2,3)
plot(y(3,1:2:127))
subplot(2,2,4)
plot(y(4,1:2:127))