24小时热门版块排行榜    

CyRhmU.jpeg
查看: 2387  |  回复: 4

nanpse

铁虫 (小有名气)

[求助] 求助 (已经-无量纲化)的非线性方程组(已经仔细检查,能算,但exitflag=0求高手

其实是后边要用matcont做常微分方程组延拓分析的。但是那个延拓软件必须有一个稳态解(充分发展的,dx/dt=0时的解)作为初值,才能开始延拓。

代码如下:已经仔细检查了一周了。能算,但是exitflag=0,估计初值太敏感。菜鸟求助啊。。。(%后边是多余的项,不用管,直接贴代码进编辑器就行。。)
小女再次感谢!

steady.m
X0=[0.278212173322273,0.912568632338739,0.085310950018308,0.912568731077048,0.087431259625095,0.087431260178122,0.499999721583934,0.000000000015842,0.087431268922953,0.053882265858466,0.087432024167792,0.087431276521867,0.500000000000000,0.087431268922956,0.087431268922956,0.500000000000000,0.087431273571452,0.500000000000001,0.078613189530066];
options=optimset('Display','iter','MaxFunEvals',100000,'MaxIter',800)
%'MaxFunEvals',100000000);
[x,fval,exitflag] =fsolve('steadys',X0,options)

steadys.m

function pp=steadys(x)
pp(1)=(1-x(1))+0.00001*48.432471*(0.0132913*x(3)*x(2)*x(7)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(5)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*1.6037e+08*2.71828^(-18.8930/x(1))*x(6)+5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(10))-2.782248*(x(1)-x(19));
pp(2)=1-x(2)-8.522*10^(-3)*1.3028e+12*2.71828^(-27.8955/x(1))*x(2)^2-5.6228*10^(-3)*1.6037e+08*2.71828^(-18.8930/x(1))*x(2)*x(3)-0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(7)-0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(6)-5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(10)-2.493604*10^(-4)*6.6735e+19*2.71828^(-45.6473/x(1))*x(9);
%-y7*c4*x(2)*x(10)
pp(3)=-x(3)+1.1722874*1.3028e+12*2.71828^(-27.8955/x(1))*x(2)^2-1.546946*1.6037e+08*2.71828^(-18.8930/x(1))*x(2)*x(3)-3.198*10^(-4)*0.01*x(3)*x(8);
%-y10*c5*x(3)*x(10);
pp(4)=1-x(4)-65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(4)-0.01923587*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(6);
pp(5)=-x(5)+2.445938*1.6037e+08*2.71828^(-18.8930/x(1))*x(2)*x(3)-5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(5)-8.496798*17.6895*2.71828^(-2.8730/x(1))*x(5)*x(8)+3.2672*10^10*1.1663e+16*2.71828^(-36.9952/x(1))*x(9)-(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(5)*x(10);
%+y16*c4*x(2)*x(10);+y18*c11
pp(6)=-x(6)-5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(6)-2.8425*10^10*1.1663e+16*2.71828^(-36.9952/x(1))*x(4)-8.496798*10^6*17.6895*2.71828^(-2.8730/x(1))*x(6)*x(8);
pp(7)=-x(7)*2.44595*10^8*1.6037e+08*2.71828^(-18.8930/x(1))*x(2)*x(4)-5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(7)-(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(7)*x(10);
%-y23*c5*x(3)*x(10);+y18*c11
pp(8)=-x(8)-5.871963*10^13*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(10)+5.05856*1.1663e+16*2.71828^(-36.9952/x(1))*x(13)+4.515847*10^5*1.1663e+16*2.71828^(-36.9952/x(1))*x(4)-134.9877*17.6895*2.71828^(-2.8730/x(1))*x(5)*x(8)+5.871963*10^13*1.1663e+16*2.71828^(-36.9952/x(1))*x(9)-8.033211*0.01*x(3)*x(8);
pp(9)=-x(9)+1.699359668*10^(-2)*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(5)-65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(9)-0.216944*6.6735e+19*2.71828^(-45.6473/x(1))*x(9);
pp(10)=-x(10)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(7)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(5)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(6)*65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(13)-x(10)*(8.496798*17.6895*2.71828^(-2.8730/x(1))*x(8)+4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1))*x(10));
%+y7*c4*x(2);+y18*c11
pp(11)=-x(11)+0.001*(0.0132913*x(3)*x(2)*x(7)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(5)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(6)+5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(10))+65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(14)-x(11)*(8.496798*17.6895*2.71828^(-2.8730/x(1))*x(8)+(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(10));
%+y7*c4*x(2);+y32*c5*x(3);+y18*c11
pp(12)=-x(12)+0.00001*(0.0132913*x(3)*x(2)*x(7)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(5)+0.0132913*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*x(6)+5.781275*1.6306e+04*2.71828^(-9.6993/x(1))*x(2)*(200*x(11)+x(10))+65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(15)-x(12)*(150.3041*2.71828^(-2.8730/x(1))*x(8)+0.036761679*2.71828^(-2.0910/x(1)))*x(10));
%+y7*c4*x(2);+y32*c5*x(3);+y18*c11
pp(13)=-x(13)*8.496798*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(10)-65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(13);
pp(14)=-x(14)+8.496798*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(11)-65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(14);
pp(15)=-x(15)+8.496798*17.6895*2.71828^(-2.8730/x(1))*x(8)*x(12)-65.3447*1.1663e+16*2.71828^(-36.9952/x(1))*x(14);
pp(16)=-x(16)*(0.5*4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(10)^2;
%+y18*c11;+x(10)*(y7*c4*x(2)+y32*c5*x(3));
pp(17)=-x(17)+(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(10)*x(11);
%+y18*c11;+x(11)*(y7*c4*x(2)+y32*c5*x(3));
pp(18)=-x(18)*(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(10)*x(12)+4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1))*x(11);
%+y18*c11;+x(12)*(y7*c4*x(2)+y32*c5*x(3))
pp(19)=(0.727148704-12.7016)*x(19)+4.897673*(x(1)-x(19));
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖


感谢参与,应助指数 +1
dbb627: 金币+1, 欢迎交流 2012-05-11 21:58:10
用1stOpt跑的结果:

x1: 0.336178337939607
x2: 1.00000000001363
x3: 1.45731439404896E-15
x4: 1
x5: -1.80058046716552E-16
x6: -7.0862384417084E-17
x7: -0.213180808520247
x8: 0.00029442958151592
x9: -1.51795582475496E-15
x10: -4.95455914027804E-12
x11: -1.94482439925578E-14
x12: 3.65096737849878E-14
x13: 0.23526045284234
x14: 4.68097286349597E-16
x15: 1.79081692185712E-14
x16: 2.7607895322144
x17: -8.50846322024266E-17
x18: -0.0647548483684659
x19: 0.0975865006697487

不知道你的公式有误问题,比如其中:
pp(16)=-x(16)*(0.5*4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x(1)))*x(10)^2;
似乎不大对吧?
2楼2012-05-11 15:34:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nanpse

铁虫 (小有名气)

引用回帖:
2楼: Originally posted by dingd at 2012-05-11 15:34:20:
用1stOpt跑的结果:

x1: 0.336178337939607
x2: 1.00000000001363
x3: 1.45731439404896E-15
x4: 1
x5: -1.80058046716552E-16
x6: -7.0862384417084E-17
x7: -0.213180808520247
x8: 0.00029442958151 ...

很谢谢呢,公式没有错呢。是阿累尼乌斯公式
还是不很对啊,都是反应的浓度,不应该有负值
还是非常谢谢呢,我回去学习一下全局的1stopt,再用maple试试
3楼2012-05-12 14:07:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★
cenwanglai: 金币+2, 谢谢回帖~ 2012-05-12 17:26:32
如果都要求大于0(开始就应该把这个条件说明):
x1: 0.336178337925227
x2: 0.999999999862475
x3: 2.27071941197048E-11
x4: 0.999999999999997
x5: 4.50305709438193E-12
x6: 1.22744048689865E-12
x7: 2.15050770012278
x8: 5.18658727927569E-14
x9: 4.62640607748565E-11
x10: 2.8267681811832E-11
x11: 2.22086150893169E-11
x12: 5.49036141029378E-12
x13: 0.695061141521546
x14: 2.9846697983009E-16
x15: 1.34476078944492E-14
x16: 0.00186288613726707
x17: 4.46468502401335E-11
x18: 0.00164117436817047
x19: 0.0975865006648326

自己带回方程验证,如果方程值均趋于0,可以说明至少在理论上是其中一组解。
4楼2012-05-12 14:17:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nanpse

铁虫 (小有名气)

引用回帖:
4楼: Originally posted by dingd at 2012-05-12 14:17:44:
如果都要求大于0(开始就应该把这个条件说明):
x1: 0.336178337925227
x2: 0.999999999862475
x3: 2.27071941197048E-11
x4: 0.999999999999997
x5: 4.50305709438193E-12
x6: 1.22744048689865E-12
x7: ...

您好!非常感谢!想请问您这个理论上能用1stopt1.5算出来吗?为什么我在1stopt里报错说“参数虽然被定义,但没有在表达式中出现呢?”

Parameter x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19;
Function(1-x1)+0.00001*48.432471*(0.0132913*x3*x2*x7+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x5+0.0132913*16306*2.71828^(-9.6993/x1)*1.6037e+08*2.71828^(-18.8930/x1)*x6+5.781275*16306*2.71828^(-9.6993/x1)*x2*x10)-2.782248*(x1-x(19))=0
1-x2-8.522*10^(-3)*1.3028*10^12*2.71828^(-27.8955/x1)*x2^2-5.6228*10^(-3)*1.6037e+08*2.71828^(-18.8930/x1)*x2*x3-0.0132913*16306*2.71828^(-9.6993/x1)*x2*x7-0.0132913*16306*2.71828^(-9.6993/x1)*x2*x6-5.781275*16306*2.71828^(-9.6993/x1)*x2*x10-2.493604*10^(-4)*6.6735e+19*2.71828^(-45.6473/x1)*x9=0
x3+1.1722874*1.3028e+12*2.71828^(-27.8955/x1)*x2^2-1.546946*1.6037e+08*2.71828^(-18.8930/x1)*x2*x3-3.198*10^(-4)*0.01*x3*x8=0
1-x4-65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x4-0.01923587*17.6895*2.71828^(-2.8730/x1)*x8*x6=0
-x5+2.445938*1.6037e+08*2.71828^(-18.8930/x1)*x2*x3-5.781275*16306*2.71828^(-9.6993/x1)*x2*x5-8.496798*17.6895*2.71828^(-2.8730/x1)*x5*x8+3.2672*10^10*1.1663e+16*2.71828^(-36.9952/x1)*x9-(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x5*x10=0
-x6-5.781275*16306*2.71828^(-9.6993/x1)*x2*x6-2.8425*10^10*1.1663e+16*2.71828^(-36.9952/x1)*x4-8.496798*10^6*17.6895*2.71828^(-2.8730/x1)*x6*x8=0
-x7*2.44595*10^8*1.6037e+08*2.71828^(-18.8930/x1)*x2*x4-5.781275*16306*2.71828^(-9.6993/x1)*x2*x7-(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x7*x10=0
-x8-5.871963*10^13*17.6895*2.71828^(-2.8730/x1)*x8*x10+5.05856*1.1663e+16*2.71828^(-36.9952/x1)*x13+4.515847*10^5*1.1663e+16*2.71828^(-36.9952/x1)*x4-134.9877*17.6895*2.71828^(-2.8730/x1)*x5*x8+5.871963*10^13*1.1663e+16*2.71828^(-36.9952/x1)*x9-8.033211*0.01*x3*x8=0
-x9+1.699359668*10^(-2)*17.6895*2.71828^(-2.8730/x1)*x8*x5-65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x9-0.216944*6.6735e+19*2.71828^(-45.6473/x1)*x9=0
-x10+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x7+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x5+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x6*65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x13-x10*(8.496798*17.6895*2.71828^(-2.8730/x1)*x8+4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1)*x10)=0
-x11+0.001*(0.0132913*x3*x2*x7+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x5+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x6+5.781275*16306*2.71828^(-9.6993/x1)*x2*x10)+65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x14-x11*(8.496798*17.6895*2.71828^(-2.8730/x1)*x8+(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x10)=0
-x12+0.00001*(0.0132913*x3*x2*x7+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x5+0.0132913*16306*2.71828^(-9.6993/x1)*x2*x6+5.781275*16306*2.71828^(-9.6993/x1)*x2*(200*x11+x10)+65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x15-x12*(150.3041*2.71828^(-2.8730/x1)*x8+0.036761679*2.71828^(-2.0910/x1))*x10)=0
-x13*8.496798*17.6895*2.71828^(-2.8730/x1)*x8*x10-65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x13=0
-x14+8.496798*17.6895*2.71828^(-2.8730/x1)*x8*x11-65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x14=0
-x15+8.496798*17.6895*2.71828^(-2.8730/x1)*x8*x12-65.3447*1.1663e+16*2.71828^(-36.9952/x1)*x14=0
-x16*(0.5*4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x10^2=0
-x17+(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x10*x11=0
-x18*(4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1))*x10*x12+4.54117*10^(-3)*8.0952*2.71828^(-2.0910/x1)*x11=0
(0.727148704-12.7016)*x(19)+4.897673*(x1-x19)=0
5楼2012-05-18 17:06:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 nanpse 的主题更新
信息提示
请填处理意见