声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 909|回复: 1

[编程技巧] 一个解微分方程组的程序

[复制链接]
发表于 2010-4-12 13:43 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
%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 编辑 ]
回复
分享到:

使用道具 举报

发表于 2010-4-14 10:22 | 显示全部楼层
LZ没有定义任何函数。。。请学习一下Matlab的函数文件的构建规则。。

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 18:32 , Processed in 0.063315 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表