24小时热门版块排行榜    

查看: 649  |  回复: 6
当前主题已经存档。

慢三儿

木虫 (小有名气)

[交流] 【求助】这个解PR方程的程序错哪里了?【已完成】

程序如下:
Tc=304.19;Pc=7.381*10^6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*10^6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=R*T/P;
V=fzero(@PReq,V0,P,T,a,b,R)


function f=PReq(V,P,T,a,b,R)
f=P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));

运行时老提示
Error in ==> SPR at 10
V=fzero(@PReq,V0,P,T,a,b,R);

why???

[ Last edited by nono2009 on 2009-10-5 at 08:11 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nono2009

超级版主 (文学泰斗)

No gains, no pains.

优秀区长优秀区长优秀区长优秀区长优秀版主

试试这个

★ ★ ★ ★ ★ ★ ★
kuhailangyu(金币+2,VIP+0):感谢double no的积极参与,o(∩_∩)o... 9-28 22:12
慢三儿(金币+5,VIP+0):谢谢!偶还想问问偶那段程序哪里错了呢?偶刚学,小菜一个...望高手指教(*^__^*) 9-29 09:14
Tc=304.19;Pc=7.381e+6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*1e6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=R*T/P;
f=@(V)P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));
V=fsolve(f,V0);
V
2楼2009-09-28 20:26:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nono2009

超级版主 (文学泰斗)

No gains, no pains.

优秀区长优秀区长优秀区长优秀区长优秀版主


kuhailangyu(金币+1,VIP+0):欢迎补充 9-28 22:12
补充一点:fzero对初值比较敏感,不推荐。
3楼2009-09-28 22:00:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

慢三儿

木虫 (小有名气)

nono2009(金币+0,VIP+0):可能还是计算的初值选得有问题。将计算条件贴出来看看。 9-29 10:23
还有 为啥分别取液态和气态的初值,计算的结果是一样滴??
4楼2009-09-29 09:49:00
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

慢三儿

木虫 (小有名气)

引用回帖:
Originally posted by 慢三儿 at 2009-9-29 09:49:
还有 为啥分别取液态和气态的初值,计算的结果是一样滴??

气态初值:V0=R*T/P
液态初值:    V0=2*b

难道是超临界,气液不分了,就出现这样的情况?
5楼2009-09-29 10:29:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nono2009

超级版主 (文学泰斗)

No gains, no pains.

优秀区长优秀区长优秀区长优秀区长优秀版主

★ ★ ★ ★ ★
慢三儿(金币+5,VIP+0):华丽丽的金币全送给no版主吧,呵呵~~版主能给诉我我那段程序错误的地方吗? 9-29 16:34
压力和温度给了,自然只有一个比容值(两相区除外)。

对超临界,一般只称超临界流体,没有气液之分。
6楼2009-09-29 11:44:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nono2009

超级版主 (文学泰斗)

No gains, no pains.

优秀区长优秀区长优秀区长优秀区长优秀版主


小木虫(金币+0.5):给个红包,谢谢回帖交流
估计是fzero的使用问题,再仔细读一下Matlab online help吧。

另外,我对比了fsolve和fzero的计算结果(如下),说明fsolve比fzero对初值的稳定性较高。

>> Tc=304.19;Pc=7.381e+6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*1e6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=R*T/P;
f=@(V)P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));
V=fzero(f,V0);
V

V =

  2.8961e-004 结果正确!

>> Tc=304.19;Pc=7.381e+6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*1e6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=R*T/P;
f=@(V)P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));
V=fsolve(f,V0);
V
Optimization terminated: norm of relative change in X is less
than max(options.TolX^2,eps) and  sum-of-squares of function
values is less than sqrt(options.TolFun).

V =

  2.8961e-004 结果正确!

>> Tc=304.19;Pc=7.381e+6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*1e6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=2*b;
f=@(V)P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));
V=fzero(f,V0);
V

V =

  2.6658e-005 结果错误!

>> f(V)

ans =

-1.1779e+019 残差很大!

>> Tc=304.19;Pc=7.381e+6;w=0.225;R=8.3145;
T=140+273.15;
P=1440/145*1e6;
Tr=T/Tc;
alfa=(1+(0.37464+1.54226*w-0.26992*w^2)*(1-Tr^0.5))^2;
ac=0.457235*(R*Tc)^2/Pc;
b=0.077796*(R*Tc)/Pc;
a=ac*alfa;
V0=2*b;
f=@(V)P-R*T/(V-b)+a/(V*(V+b)+b*(V-b));
V=fsolve(f,V0);
V
Optimization terminated: norm of relative change in X is less
than max(options.TolX^2,eps) and  sum-of-squares of function
values is less than sqrt(options.TolFun).

V =

  2.8961e-004 结果正确!

[ Last edited by nono2009 on 2009-9-30 at 00:23 ]
7楼2009-09-30 00:21:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 慢三儿 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] FileCode能看出啥? +7 要乐观耀哥 2026-08-10 16/800 2026-08-11 08:18 by wmfsnow
[基金申请] 小木虫上这么多卖论文的,真有人买论文么?感觉没必要啊 +4 Tide man 2026-08-10 4/200 2026-08-11 07:45 by emeraldine
[基金申请] 基金中了 +15 laoda193707 2026-08-06 15/750 2026-08-11 00:11 by jiafei2190
[基金申请] 确定了,国自然21号放榜 +6 布布和一二 2026-08-10 7/350 2026-08-10 19:15 by 2000zf36392
[基金申请] 关于代码变化问题,想知道的进来 +17 且听虎啸 2026-08-07 24/1200 2026-08-10 18:35 by zhangduo2008
[基金申请] 关于Filecode分析方法 +3 majunge000 2026-08-10 3/150 2026-08-10 15:46 by lch2012
[基金申请] 奇怪,两个人的filecode固定段从头到尾一模一样 +4 布布和一二 2026-08-10 4/200 2026-08-10 14:50 by 医学老男孩
[基金申请] 好奇怪的filecode +4 布布和一二 2026-08-08 5/250 2026-08-10 14:20 by 冰心玉壶晴
[基金申请] 2026国自然放榜时间 +9 布布和一二 2026-08-08 9/450 2026-08-10 11:22 by xxxx2020
[基金申请] 关于filecode,很负责任的告诉大家 +6 爱看书的可乐 2026-08-08 7/350 2026-08-08 22:13 by a_niu
[基金申请] 关于豆爷回答的JTJC与%2F数量 +5 yang182083 2026-08-06 7/350 2026-08-08 18:29 by zhanghaozhu
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 化学口download_prp&fileCode的固定段好像这几天一直没变,有变的大神么? +3 Tide man 2026-08-07 4/200 2026-08-07 22:39 by Tide man
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] 听说今天filecode变了 +24 布布和一二 2026-08-06 47/2350 2026-08-07 16:02 by zhiyanjiang
[基金申请] filecode变化情况 +6 布布和一二 2026-08-07 22/1100 2026-08-07 14:45 by 且听虎啸
[基金申请] filecode +14 等待解的谜 2026-08-06 19/950 2026-08-07 12:20 by wlwhappy
[基金申请] filecode +8 布布和一二 2026-08-06 11/550 2026-08-06 20:41 by tangpu318
[基金申请] 求各位大神看下 100+6 hpkpkpkp 2026-08-05 33/1650 2026-08-06 14:49 by zhiyanjiang
[基金申请] 好消息?这个有何含义??? +8 Tide man 2026-08-05 10/500 2026-08-05 16:14 by xmuxiaoyu
信息提示
请填处理意见