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

[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开始）。