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

[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