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

[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);

[yuyiyy10103]

补充一下，后面的那张图片就是那１８个公式的方程式，看起来可能会方便理解一点。

还有，我用runge-kutta四阶法求近似，

这里是runge-kutta的程序: rk4.m

function [x, y, dy] = rk4(deriv,n,x,dx,y)

x0 = x;
y0 = y;

[y,dy1] = feval(deriv,x0,y); % k1=f(x0,y0)
for i = 1:n
y(i) = y0(i) + 0.5*dx*dy1(i);
end

xm = x0 + 0.5*dx; % x0+0.5h

[y,dy2] = feval(deriv,xm,y);
% k2=f(x0+0.5h,y0+0.5k1*h)
for i = 1:n
y(i) = y0(i) + 0.5*dx*dy2(i); % x0+0.5*h, y0+0.5k2*h
end
[y,dy3] = feval(deriv,xm,y);
% k3=f(x0+0.5*h, y0+0.5k2*h)
for i = 1:n
y(i) = y0(i) + dx*dy3(i); % x0+h, y0+0.5k3*h
end
x = x0 + dx;
[y,dy] = feval(deriv,x,y);
% k4=f(x0+h,y0+0.5k3*h)
for i = 1:n
dy(i) = (dy1(i) + 2*(dy2(i) + dy3(i)) + dy(i))/6;
y(i) = y0(i) + dx*dy(i); % y
end

[yuyiyy10103]

这是第一步求气缸体积vc,ve的函数:
vclc = 8.000e-06;
vcle = 1.000e-05;
b1 = 3.540e-02;
b2 = 3.540e-02;
crank = 8.500e-03;
dcomp = 5.590e-02;
dexp = 5.590e-02;
% post traitement (geometrie)
acomp = pi*dcomp^2/4.0;
aexp = pi*dexp^2/4.0;
yoke = sqrt(b1^2 + b2^2);
ymax = sqrt((yoke + crank)^2 - b2^2);
ymin = sqrt((yoke - crank)^2 - b2^2);
vswc = acomp*(ymax - ymin); % compression swept volumes [m^3]
vswe = aexp*(ymax - ymin); % expansion swept volumes [m^3]
% partie de fonction volume.m
theta = [0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360];
for i=1:37
sinth(i) = sin(theta(i));
costh(i) = cos(theta(i));
bth(i) = (b1^2 - (crank*costh(i))^2)^0.5;
ye(i) = crank*(sinth(i) + (b2/b1)*costh(i)) + bth(i);
yc(i) = crank*(sinth(i) - (b2/b1)*costh(i)) + bth(i);
ve(i) = vcle + aexp*(ye(i) - ymin);
vc(i) = vclc + acomp*(yc(i) - ymin);
dvc(i) = acomp*crank*(costh(i) + (b2/b1)*sinth(i) + crank*sinth(i)*costh(i)/bth(i));
dve(i) = aexp*crank*(costh(i) - (b2/b1)*sinth(i) + crank*sinth(i)*costh(i)/bth(i));
fprintf('dve= %.8f\n',dve);
fprintf('dvc= %.8f\n',dvc);
fprintf('vc= %.8f\n',vc);
fprintf('ve= %.8f\n',ve);
end
plot(theta,vc)

[yuyiyy10103]

付给初始值的函数: adiab.m

epsilon = 1; % Allowable error in temerature (K)
max_iteration = 20; % Maximum number of iterations to convergence
step = 10; % for saving values in var, dvar matrices
dtheta = 1; % integration increment (radians)
% Initial conditions:
y(the) = th;
y(tck) = tk;
y(te) = th;
y(tc) = tk;
iter = 0;
terror = 10*epsilon; % Initial error to enter the loop
% Iteration loop to cyclic convergence
while ((terror >= epsilon)&(iter < max_iteration))
% cyclic initial conditions
tc0 = y(tc);
te0 = y(te);
theta = 0;
y(qk) = 0;
y(qr) = 0;
y(qh) = 0;
y(wc) = 0;
y(we) = 0;
y(w) = 0;
fprintf('iteration %d: tc = %.1f[K], te = %.1f[K]\n',iter,y(tc),y(te))
for(i = 1:1:360)
[theta,y,dy] = rk4('dadiab',7,theta,dtheta,y);
end
terror = abs(tc0 - y(tc)) + abs(te0 - y(te));
iter = iter + 1;
end

if (iter >= max_iteration)
fprintf('No convergence within %d iteration\n',max_iteration)
end
for(j = 1:1:10)
[theta,y,dy] = rk4('dadiab',7,theta,dtheta,y);
end

[yuyiyy10103]

主程序最后一段我调用adiab.m dadiab.m rk4.m

老说我没有定义vc什么的!!!!

其实就是调用rk4.m的过程中总是出问题，不知道怎么调用。。。

高手帮帮忙，谢谢了!

[watcher]

如果按
rk4('dadiab',7,theta,dtheta,y);
这种格式调用dadiab，那么dadiab好像应该是函数，即以function开头

lz的做了很多工作，但是遗憾得告诉你，matlab中ode系列函数可以完善的解决你的问题。也就是说，除非你是故意想锻炼一下你的matlab编程技巧，完全没必要自己编写runge-kutta求解函数。
建议先看一下matlab中关于ode求解的帮助文档（通常我会从函数ode45开始）。