超难解微分方程组的一个程序

[yuyiyy10103]

总共有１８个方程式，

最后需要解出tc,te

现在已知的有vk,vh,vr,th,tk,tr是固定的，

vc,ve 是两个随轴承角度theta而变化的气缸体积的矩阵，

这是 方程列表的 函数: dadiab.m

vot = vk/tk + vr/tr + vh/th;
y(p) = (mgas*rgas/(y(vc)/y(tc) + vot + y(ve)/y(te)));
top = -y(p)*(dy(vc)/y(tck) + dy(ve)/y(the));
bottom = (y(vc)/(y(tck)*gama) + vot + y(ve)/(y(the)*gama));
dy(p) = top/bottom;

% mass accumulations and derivatives:
y(mc) = y(p)*y(vc)/(rgas*y(tc));
y(mk) = y(p)*vk/(rgas*tk);
y(mr) = y(p)*vr/(rgas*tr);
y(mh) = y(p)*vh/(rgas*th);
y(me) = y(p)*y(ve)/(rgas*y(te));
dy(mc) = (y(p)*dy(vc) + y(vc)*dy(p)/gama)/(rgas*y(tck));
dy(me) = (y(p)*dy(ve) + y(ve)*dy(p)/gama)/(rgas*y(the));
dpop = dy(p)/y(p);
dy(mk) = y(mk)*dpop;
dy(mr) = y(mr)*dpop;
dy(mh) = y(mh)*dpop;

% mass flow between cells:
y(mck) = -dy(mc);
y(mkr) = y(mck) - dy(mk);
y(mhe) = dy(me);
y(mrh) = y(mhe) + dy(mh);

% conditional temperatures between cells:
y(tck) = tk;
if((ymck)>0)
y(tck) = y(tc);
end
y(the) = y(te);
if(y(mhe)>0)
y(the) = th;
end

% 7 derivatives to be integrated by rk4:
% working space temperatures:
dy(tc) = y(tc)*(dpop + dy(vc)/y(vc) - dy(mc)/y(mc));
dy(te) = y(te)*(dpop + dy(ve)/y(ve) - dy(me)/y(me));

% energy:
dy(qk) = vk*dy(p)*cv/rgas - cp*(y(tck)*y(mck) - tk*y(mkr));
dy(qr) = vr*dy(p)*cv/rgas - cp*(tk*y(mkr) - th*y(mrh));
dy(qh) = vh*dy(p)*cv/rgas - cp*(th*y(mrh) - y(the)*y(mhe));
dy(wc) = y(p)*dy(vc);
dy(we) = y(p)*dy(ve);

% Net work done:
dy(w) = dy(wc) + dy(we);
y(w) = y(wc) + y(we);