应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 怎样设置求解精度，一直提醒我方程解接近于零 应该怎么改呢
51/1返回列表
查看: 683  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

杜�?Ramos

铁虫 (小有名气)

[交流] 怎样设置求解精度，一直提醒我方程解接近于零 应该怎么改呢 已有2人参与

function equation()
global Rg M T Veff P
Rg=8.31;
M=5.8*10^-3;
T=273.15;
Veff=1;
P=10000;
n0=[53.156*10^-3 0 0 5.98*10^-9 0 0 0 0];
n=fsolve(@myfun,n0);
end

function F=myfun(n)
global Rg M T Veff P
n1=n(1);
n2=n(2);
n3=n(3);
n4=n(4);
n5=n(5);
n6=n(6);
n7=n(7);
n8=n(8);
kf1=1.912*10^(23)*(T^(-1.83))*exp(-1.185*10^5/(8.31*T));
kr1=2.2*10^(22)*(T^(-2))*exp(0/(8.31*T));
kf2=4.515*10^(17)*(T^(-0.64))*exp(-1.189*10^5/(8.31*T));
kr2=6.165*10^(15)*(T^(-0.5))*exp(0/(8.31*T));
kf3=2.97*10^(6)*(T^(2.02))*exp(-1.34*10^4/(8.31*T));
kr3=1.465*10^(5)*(T^(2.11))*exp(2.904*10^3/(8.31*T));
kf4=1.475*10^(12)*(T^(0.6))*exp(0/(8.31*T));
kr4=3.09*10^(12)*(T^(0.53))*exp(-4.887*10^4/(8.31*T));
kf5=2.667*10^(4)*(T^(2.65))*exp(-4.88*10^3/(8.31*T));
kr5=3.82*10^(12)*(T^(0))*exp(-7.948*10^3/(8.31*T));
kf6=2.298*10^(9)*(T^(1.4))*exp(-1.832*10^4/(8.31*T));
kr6=2.16*10^(8)*(T^(1.52))*exp(-3.45*10^3/(8.31*T));
kf7=3.25*10^(13)*(T^(0))*exp(0/(8.31*T));
kr7=3.252*10^(12)*(T^(0.33))*exp(-5.328*10^4/(8.31*T));
kf8=7.079*10^(13)*(T^(0))*exp(-2.95*10^2/(8.31*T));
kr8=2.027*10^(10)*(T^(-0.72))*exp(-3.684*10^4/(8.31*T));
kf9=1.838*10^(10)*(T^(0.59))*exp(-3.089*10^4/(8.31*T));
kr9=1*10^(12)*(T^(0))*exp(0/(8.31*T));
kf10=8.66*10^(3)*(T^(2.68))*exp(-1.856*10^4/(8.31*T));
kr10=9.550*10^(6)*(T^(2))*exp(-3.97*10^3/(8.31*T));
kf11=1.66*10^(13)*(T^(0))*exp(-8.23*10^2/(8.31*T));
kr11=3.164*10^(12)*(T^(0.35))*exp(-5.551*10^4/(8.31*T));
kf12=1*10^(14)*(T^(-0.37))*exp(0/(8.31*T));
kr12=2.951*10^(14)*(T^(0))*exp(-4.843*10^4/(8.31*T));
kf13=4.2*10^(14)*(T^(0))*exp(-1.198*10^4/(8.31*T));
kr13=4.624*10^(16)*(T^(-0.35))*exp(-5.067*10^4/(8.31*T));
nt=n1+n2+n3+n4+n5+n6+n7+n8;
eq1=-(nt*Rg*T)/P*((kf1*(n1/Veff)*(M/Veff)-kr1*(n2/Veff)*(n3/Veff)*(M/Veff))+(kf3*(n1/Veff)*(n5/Veff)-kr3*(n3/Veff)*(n3/Veff))+(kf6*(n1/Veff)*(n2/Veff)-kr6*(n7/Veff)*(n3/Veff))+(kf9*(n1/Veff)*(n6/Veff)-kr9*(n8/Veff)*(n3/Veff)));
eq2=-(nt*Rg*T)/P*((kr1*(n2/Veff)*(n3/Veff)*(M/Veff)-kf1*(n1/Veff)*(M/Veff))+(kf4*(n2/Veff)*(n4/Veff)*(M/Veff)-kr4*(n6/Veff)*(M/Veff))+(kf5*(n2/Veff)*(n3/Veff)-kr5*(n5/Veff)*(n7/Veff))+(kf6*(n2/Veff)*(n1/Veff)-kr6*(n3/Veff)*(n7/Veff))+(kf8*(n2/Veff)*(n6/Veff)-kr8*(n3/Veff)*(n3/Veff))+(kf11*(n6/Veff)*(n2/Veff)-kr11*(n7/Veff)*(n4/Veff)));
eq3=-(nt*Rg*T)/P*((kr1*(n2/Veff)*(n3/Veff)*(M/Veff)-kf1*(n1/Veff)*(M/Veff))+(2*kr3*(n3/Veff)*(n3/Veff)-kr3*(n1/Veff)*(n5/Veff))+(kf5*(n2/Veff)*(n3/Veff)-kr5*(n5/Veff)*(n7/Veff))+(kr6*(n3/Veff)*(n7/Veff)-kf6*(n2/Veff)*(n1/Veff))+(kr7*(n3/Veff)*(n4/Veff)-kf7*(n6/Veff)*(n5/Veff))+(2*kr8*(n3/Veff)*(n3/Veff)-kf8*(n6/Veff)*(n2/Veff))+(kr9*(n8/Veff)*(n3/Veff)-kf9*(n1/Veff)*(n6/Veff))+(kf10*(n3/Veff)*(n6/Veff)-kr10*(n8/Veff)*(n5/Veff))+(2*kf12*(n3/Veff)*(n3/Veff)*(M/Veff)-kr11*(n8/Veff)*(M/Veff)));
eq4=-(nt*Rg*T)/P*((kf2*(n4/Veff)*(M/Veff)-kr2*(n5/Veff)*(n5/Veff)*(M/Veff))+(kf4*(n2/Veff)*(n4/Veff)*(M/Veff)-kr4*(n6/Veff)*(M/Veff))+(kr7*(n3/Veff)*(n4/Veff)-kf7*(n6/Veff)*(n5/Veff))+(kr11*(n7/Veff)*(n4/Veff)-kf11*(n6/Veff)*(n2/Veff))+(kr13*(n8/Veff)*(n4/Veff)-kf13*(n6/Veff)*(n6/Veff)));
eq5=-(nt*Rg*T)/P*((2*kr2*(n5/Veff)*(n5/Veff)*(M/Veff)-kf2*(n4/Veff)*(M/Veff))+(kf3*(n1/Veff)*(n5/Veff)-kr3*(n3/Veff)*(n3/Veff))+(kr5*(n5/Veff)*(n7/Veff)-kf5*(n2/Veff)*(n3/Veff))+(kf7*(n6/Veff)*(n5/Veff)-kr7*(n3/Veff)*(n4/Veff))+(kr10*(n8/Veff)*(n5/Veff)-kf10*(n3/Veff)*(n6/Veff)));
eq6=-(nt*Rg*T)/P*((kr4*(n6/Veff)*(M/Veff)-kf4*(n2/Veff)*(n4/Veff)*(M/Veff))+(kf7*(n6/Veff)*(n5/Veff)-kr7*(n3/Veff)*(n4/Veff))+(kf8*(n6/Veff)*(n2/Veff)-kr8*(n3/Veff)*(n3/Veff))+(kf9*(n1/Veff)*(n6/Veff)-kr9*(n3/Veff)*(n8/Veff))+(kf10*(n3/Veff)*(n6/Veff)-kr10*(n8/Veff)*(n5/Veff))+(kf11*(n6/Veff)*(n2/Veff)-kr11*(n7/Veff)*(n4/Veff))+(2*kf13*(n6/Veff)*(n6/Veff)-kr13*(n8/Veff)*(n4/Veff)));
eq7=-(nt*Rg*T)/P*((kr5*(n5/Veff)*(n7/Veff)-kf5*(n2/Veff)*(n3/Veff))+(kr11*(n4/Veff)*(n7/Veff)-kf11*(n6/Veff)*(n2/Veff)));
eq8=-(nt*Rg*T)/P*((kr9*(n3/Veff)*(n8/Veff)-kf9*(n1/Veff)*(n6/Veff))+(kr10*(n8/Veff)*(n5/Veff)-kf10*(n3/Veff)*(n6/Veff))+(kr12*(n8/Veff)*(M/Veff)-kf12*(n3/Veff)*(n3/Veff)*(M/Veff))+(kr13*(n8/Veff)*(n4/Veff)-kf13*(n6/Veff)*(n6/Veff)));
F=[eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8];
end
回复此楼

» 猜你喜欢
1楼 2017-10-23 14:21:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

moguowei00

铜虫 (小有名气)
好难

发自小木虫Android客户端
回复此楼
2楼2017-10-23 14:34:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

杜�?Ramos

铁虫 (小有名气)
引用回帖:
2楼: Originally posted by moguowei00 at 2017-10-23 14:34:08
好难

函数表达式有点长，就是求解属于很简单。我不知道怎样设置求解的取值精度
赞一下 回复此楼
3楼2017-10-23 14:47:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

pdl9527

专家顾问 (小有名气)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
你的方程的形式是 ax-b=0， 所以这里提示你最后的结果能够保证这个差值接近0了，如果你把求解的那一行改成 [n,fval,exitflag,output]=fsolve(@myfun,n0)， 你会发现exitflag=3，一般来说exitflag>0的情况下求得的解都是可靠的。exitflag具体代表的意思你可以在命令框输入doc fsolve自行了解。
赞一下 回复此楼
4楼2017-10-23 16:52:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

杜�?Ramos

铁虫 (小有名气)
引用回帖:
4楼: Originally posted by pdl9527 at 2017-10-23 16:52:40
你的方程的形式是 ax-b=0， 所以这里提示你最后的结果能够保证这个差值接近0了，如果你把求解的那一行改成 =fsolve(@myfun,n0)， 你会发现exitflag=3，一般来说exitflag>0的情况下求得的解都是可靠的。exitf ...

谢谢 我试一下
回复此楼
5楼2017-10-23 17:11:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 杜�?Ramos 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭