我正在尝试求解 6 个 ODE 方程和 1 个 DAE 方程。ODE方程已在空间域进行了离散化，并使用了ode15s MATLAB求解器对时域方程进行了求解。我已经定义了一个对角质量矩阵，包括 ODE 的单位元素和 DAE 方程的零值。假设“n”是沿空间的节点数，M 定义如下：

  O=ones(1,6*n);
  Ze=zeros(1,n);
  dia=horzcat(O,Ze);
  M=diag(dia);

这是我的主要和功能 m.files 的重要部分：

.
.
.

y0(1:n)=zeros(n,1); %C0
y0(n+1:2*n)=zeros(n,1); %qc0
y0(2*n+1:3*n)=Tini*ones(n,1); %T0
y0(3*n+1:4*n)=Tini*ones(n,1); %Ts0
y0(4*n+1:5*n)=zeros(n,1); %N0
y0(5*n+1:6*n)=zeros(n,1); %qN0
y0(6*n+1:7*n)=zeros(n,1); %U0
t_int=1;
tf=3000;
nt=(tf/t_int);
tspan=(1:t_int:tf);
options=odeset('Mass',M);
[t,y]=ode15s(@adsor,tspan,y0,options);

和

function dydt=adsor(t,y)
.
.
.
dydt=zeros(7*n,1);
y(1,1)=Cin(1); %BC of C
y(n,1)=y(n-1,1); %BC of C
y(2*n+1,1)=Tg; %BC of Tg
y(3*n,1)=y(3*n-1,1); %BC of Tg
y(4*n+1,1)=Cin(2); %BC of N
y(5*n,1)=y(5*n-1,1); %BC of N
y(6*n+1,1)=Q/A; %BC of U
C(1:n,1)=y(1:n,1);
qc(1:n,1)=y(n+1:2*n,1);
T(1:n,1)=y(2*n+1:3*n,1);
Ts(1:n,1)=y(3*n+1:4*n,1);
N(1:n,1)=y(4*n+1:5*n,1);
qN(1:n,1)=y(5*n+1:6*n,1);
U(1:n,1)=y(6*n+1:7*n,1);

for i=2:n-1

    dqcdt(i,1)=K1*(qsc-qc(i));
    dCdt(i,1)=aEC(i).*C(i+1)+aWC(i).*C(i-1)-aPC(i).*C(i)+dqcdt(i-1));
    dqNdt(i,1)=K2*(qcn-qN(i));
    dNdt(i,1)=aEC(i).*N(i+1)+aWC(i).*N(i-1)-aPC(i).*N(i))+dqNdt(i-1);
    dTdt(i,1)=aET(i).*T(i+1)+aWT(i).*T(i-1)-aPT(i)+*T(i))+Ts(i);
    dTsdt(i,1)=T(i)-Ts(i)+dqcdt(i)+dqNdt(i);
    dUdt(i,1)=U(i)-U(i-1)+dqcdt(i)+dqNdt(i);

end
dydt(1:7*n,1)=[dCdt;0 ;dqcdt;0 ;dTdt;0 ;dTsdt;0 ;dNdt;0 ;dqNdt;0 ;dUdt;0];

位于最后一行的方程（即 dUdt(i,1)=...）是由定义的质量矩阵识别的 DAE 方程（请参见前面定义的“M”参数）。我不确定这种方法是否正确指定 DAE 方程！无论如何，当我运行这些代码时，我面临以下错误：

使用 daeic12 时出错（第 77 行）此 DAE 的索引似乎大于 1。

我检查了雅可比矩阵，发现它是 t0 处的*n 矩阵。有两列（第一列和最后一列）都为零，我认为它使雅可比矩阵奇异（如https://ch.mathworks.com/help/matlab/math/solve-differential-algebraic-方程-daes.html）。

您能否指导我如何克服此错误并修复我的代码？任何帮助表示赞赏。

最好的祝福