一个解微分方程组的程序

cjw718

%main.f
global t1 p1 t2 p2 p3 q1 tj vs mj cf a b hh v n h2 y0
clear
t1=244;
p1=10.9;
t2=293;
p2=10.8;
P3=10.7;
tj=402;
q1=68000;
vs=210;mj=1603900;cj=0.46;cf=1.49;
a=0;b=300;
hh=0.1;
alfa=0.1;
[p,t,h2,s,v,x,r] = WASPCHS('PT',p2,t2);
n=3000;
y0=[1./v;h2;p2;tj];
fid=fopen('e:\matlab\shuju.txt','wt')
for i=1:n
    if (i>=0)&(i<=2400)
        q1=q1*(1+alfa/600);
    elseif (i>2400)&&(i<=3000)
        continue;
    end   
[x,y]=RungeKutta4('func',a,b,y0,n)
   
fprintf(fid,'5%10.4f',th_1,y(1),y(2),y(3),y(4));
end
fclose(fid);

%RungeKutta4.m
global t1 p1 t2 p2 p3 q1 tj vs mj cf a b hh v h2 n y0
hh=(b-a)/n;

x(1)=0;
for ii=1:n
  x(ii+1)=x(ii)+hh;
  k1=feval('func',x(ii),y(:,ii));
  k2=feval('func',x(ii)+hh/2,y(:,ii)+hh*k1/2);
  k3=feval('func',x(ii)+hh/2,y(:,ii)+hh*k2/2);
  k4=feval('func',x(ii)+hh,y(:,ii)+hh*k3);
  y(:,ii+1)=y(:,ii)+hh*(k1+2*k2+2*k3+k4)/6;
end

%func.m
global t1 p1 t2 p2 p3 q1 tj vs mj cf a b hh v y0
format long;
v1=WASPCHS('PT2V',p1,t1);
ro1=1./v1;
p=p1+0.001;
v2=WASPCHS('PT2V',p,t1);
ro2=1./v2;
drodp=(ro2-ro1)./0.001;
[p,t,h1,S,v3,X,R] = WASPCHS('PT',p1,t1);
ro3=1./v3;
h=h1+1;
v4=WASPCHS('PH2V',p1,h1);
ro4=1./v4;
drodh=ro4-ro3;
din =(1./cf.*ro1.*(p1-y(3)).*1000)^0.5;
dout=(1./cf.*y(2).*(p2-p3)*1000).^0.5;
q2=6.7*din.^0.8.*(y(4)-t2);
dy=zero(4,1);
dy(1)=(din-dout)./vs;
dy(2)=((din-dout)-vs.*drodp.*dy(3))./vs./drodh;
dy(3)=((y(1)./drodp+y(2)).*(din-dout)-(din.*h1-dout.*y(2)+q2))./vs/(1+drodp./drodh.*y(1));
dy(4)=(q1-q2)./mj./cj;  

运行后显示：
??? One or more output arguments not assigned during call to 'E:\MATLAB\chengxu\RungeKutta4.m (RungeKutta4)'.
Error in ==> MF at 24
[x,y]=RungeKutta4('func',a,b,y0,n)
我估计是程序没有执行FUNC.M这里面的函数。

由于是刚学的M文件，可能会有很多错误，看了好几天了，也不知道哪儿错了，请大家帮忙改一下吧，谢谢了。

[ 本帖最后由 cjw718 于 2010-4-12 16:24 编辑 ]

maigicku

LZ没有定义任何函数。。。请学习一下Matlab的函数文件的构建规则。。