| 查看: 614 | 回复: 2 | ||
[求助]
求解非线性六元方程组，时间紧，自己来不及学了，麻烦大家帮忙，问题简单，悬赏多。 已有1人参与
|
|
用MATLAB 非线性求解，不会 >> syms qn dn dr pn pr cn cr sn sr hn hr t e a o A T u1 k Fn Fr zn zr >>o=0.8; >>a=0.1; >>cn=0.4; >>cr=0.2; >>k=2; >>T=1; >>hn=0.1; >>sn=0.1; >>hr=0.05; >>sr=0.05; >>A=0.1; >> dn=1-(pn-pr)/(1-o)>> dr=(o*pn-pr)/(o-o^2) >> Fn=20*(qn-dn) >> Fr=20*(qn*e*(a*t+1)-dr) >> zn=qn-dn >> zr=qn*e*(a*t+1)-dr >> eq1='-(pn+hn-cn-sn)*Fn+pn+hn-cn-e*(a*t+1)*(pr+hr+A*(T-t)-cn-sr)*Fr+e*(a*t+1)*(pr+hr+A*(T-t)-cn)+u1*(1-e*(a*t+1))=0' >> eq2='qn+(hn-hr)/(1-o)+(1+qn-(2*pn-pr-sn-cn+hn)/(1-o))*Fn+Fr*(pr+hr+A*(T-t)-cn-sr)/(1-o)+10*zn^2=0' >> eq3='e*qn*(a*t+1)+(hr-o*hn)/(o-o^2)+Fn*(pn+hn-cn-sn)/(1-o)+Fr*(o*pn-2*pr+cn+sr-A*(T-t)-hr-o*(1-o)*e*qn*(a*t+1))/(o-o^2)+10*0.05^2-10*zr^2=0' >> eq4='a*e*qn*(hr+pr+A*(T-2*t-1/a)-cn)+a*e*qn*(sr+cn-hr-A*(T-2*t-1/a)-pr)*Fr-A*(20*dr*zr+10*zr^2)-u1*qn*e*a=0' >> eq5='qn*(a*t+1)*(pr+hr+A*(T-t)-cn)-qn*(a*t+1)*(pr+hr+A*(T-t)-cn-sr)*Fr-2*k*e-u1*a*t*qn=0' >> eq6='qn*(1-e*(a*t+1))=0' 最终求解qn pn pr t e u1 是用fsolve算吗，初始值不知道怎么赋，大致推迟大概是0.5 0.5 0.5 0.5 0.5 谢谢 |
回复此楼» 猜你喜欢
听说今天filecode变了
已经有34人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
基金中了
已经有6人回复
-
filecode
已经有14人回复
-
固定端突然变了，今天
已经有7人回复
-
filecode
已经有11人回复
-
大家散了吧，后缀研究没有意义，别浪费时间了，过好目前的每一天，不要焦虑
已经有6人回复
-
一个有机合成实验室都需要哪些设备？
已经有12人回复
-
求各位大神看下
已经有33人回复
-
咨询面上基金
已经有7人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 求助：非线性方程组的求解（郁闷中）
已经有26人回复
- 数学小白求助，解超定非线性方程组，求大神帮忙
已经有4人回复
- Matlab用牛顿法求解非线性方程组问题
已经有8人回复
- 用matlab求解非线性方程组说无解，一定是方程组本身无解，还是有可能程序有问题呢？
已经有11人回复
- matlab如何求解一个非线性微分方程组
已经有8人回复
- MATLAB求解非线性方程组
已经有5人回复
- 关于matlab求解非线性指数方程组出现问题
已经有3人回复
- matlab求解非线性方程组，错误提示怎么解决
已经有5人回复
- matlab求解非线性方程组，求助！
已经有6人回复
- 求助！matlab用fsolve函数求解非线性方程组的问题！
已经有19人回复
- 如何使用matlab求解非线性方程组的所有整数解？
已经有9人回复
- matlab求解非线性方程组
已经有16人回复
- 求高人指点用matlab求解非线性方程组，解决了追加100金币；
已经有11人回复
- matlab的fsove 命令求解非线性方程组
已经有6人回复
- 【求助】用mathematica 5.0求解一个非线性方程组失败，特发帖求助！
已经有5人回复
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
旖旎落下: 金币+50, ★★★很有帮助, 好。 2014-01-09 10:23:14
感谢参与，应助指数 +1
旖旎落下: 金币+50, ★★★很有帮助, 好。 2014-01-09 10:23:14
|
以下是如何求数值解。注意，为了简化，让优化工具箱自行计算方程组的数值一阶微分（雅可比矩阵）和二阶微分（Hessian矩阵）。如果结果不理想，可以先调节 options 中的参数。如果还不理想，可以利用符号运算工具箱求出雅可比矩阵和 Hessian 矩阵。 1. 建立一个函数文件 func1.m function y=func1(x) %未知量 qn=x(1); pn=x(2); pr=x(3); t=x(4); e=x(5); u1=x(6); %常量 o=0.8; a=0.1; cn=0.4; cr=0.2; k=2; T=1; hn=0.1; sn=0.1; hr=0.05; sr=0.05; A=0.1; %简化表达式 dn=1-(pn-pr)/(1-o); dr=(o*pn-pr)/(o-o^2); Fn=20*(qn-dn); Fr=20*(qn*e*(a*t+1)-dr); zn=qn-dn; zr=qn*e*(a*t+1)-dr; %非线性方程组 y = [-(pn+hn-cn-sn)*Fn+pn+hn-cn-e*(a*t+1)*(pr+hr+A*(T-t)-cn-sr)*Fr+e*(a*t+1)*(pr+hr+A*(T-t)-cn)+u1*(1-e*(a*t+1)); qn+(hn-hr)/(1-o)+(1+qn-(2*pn-pr-sn-cn+hn)/(1-o))*Fn+Fr*(pr+hr+A*(T-t)-cn-sr)/(1-o)+10*zn^2; e*qn*(a*t+1)+(hr-o*hn)/(o-o^2)+Fn*(pn+hn-cn-sn)/(1-o)+Fr*(o*pn-2*pr+cn+sr-A*(T-t)-hr-o*(1-o)*e*qn*(a*t+1))/(o-o^2)+10*0.05^2-10*zr^2; a*e*qn*(hr+pr+A*(T-2*t-1/a)-cn)+a*e*qn*(sr+cn-hr-A*(T-2*t-1/a)-pr)*Fr-A*(20*dr*zr+10*zr^2)-u1*qn*e*a; qn*(a*t+1)*(pr+hr+A*(T-t)-cn)-qn*(a*t+1)*(pr+hr+A*(T-t)-cn-sr)*Fr-2*k*e-u1*a*t*qn; qn*(1-e*(a*t+1))]; 2. 建立一个主程序 solve_6_unknowns.m %清屏，清工作区 clc clear all %设置优化算法参数 maxiter = 20000; maxfuneval = length(x0)*maxiter; options = optimset('Display‘， ’off',... 'GradObj', 'off',... 'Hessian', 'off',... 'TolX', 1e-6,... 'TolFun', 1e-6,... 'MaxIter', maxiter,... 'MaxFunEvals', maxfuneval); %设置初始值 x0 = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5]; %调用优化函数 [x, fval, exitflag, output] = fsolve(@func1, x0, options); %取结果 qn=x(1); pn=x(2); pr=x(3); t=x(4); e=x(5); u1=x(6); %显示结果 disp(['Iterations: ',num2str(output.iterations)]) disp(['Func Evals: ', num2str(output.funcCount)]) disp(['Algorithm: ',output.algorithm]) disp(['exit flag = ',num2str(exitflag)]) disp(['error = ( ',num2str(fval','%-15.6e'),' )']) disp(['qn = ( ',num2str(qn,'%-15.6f'),' )']) disp(['pn = ( ',num2str(pn,'%-15.6f'),' )']) disp(['pr = ( ',num2str(pr,'%-15.6f'),' )']) disp(['t = ( ',num2str(t,'%-15.6f'),' )']) disp(['e = ( ',num2str(e,'%-15.6f'),' )']) disp(['u1 = ( ',num2str(u1,'%-15.6f'),' )']) |
2楼2014-01-08 22:57:57
3楼2014-01-09 10:23:46
50