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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 用1stOpt解出非线性方程组的解作为初值用Matlab求解，解不出来
51/1返回列表
查看: 1623  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

mjjackey

新虫 (初入文坛)

[求助] 用1stOpt解出非线性方程组的解作为初值用Matlab求解，解不出来 已有1人参与

用1stOpt解出非线性方程组的解作为初值用Matlab求解，解不出来
这是一个非线性方程组（Matlab表示）：
CODE:
function  f= solveMagneticDipoleEquation(x)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
%已知数
a=-150;
b=150;
Bx=-48030;
By=29370;
Bz=-19077.3;
%未知数
Bt=x(1);
c=x(2);
m=x(3);
n=x(4);
p=x(5);
Rl=sqrt(a.^2+b.^2+c.^2);
f(1)=Bt*(3*(m*a+n*b+p*c)*a/Rl.^5-m/Rl.^3)/Bx-1;
f(2)=Bt*(3*(m*a+n*b+p*c)*b/Rl.^5-n/Rl.^3)/By-1;
f(3)=Bt*(3*(m*a+n*b+p*c)*c/Rl.^5-p/Rl.^3)/Bz-1;
f(4)=m.^2+n.^2+p.^2-1;
end

用一开始建模型的人给的粗略的初始值范围，用默认设置选项，解出的解明显不对，警告提示要用Levenberg-Marquardt算法，后改用这种算法options=optimset('Display','iter','Algorithm','Levenberg-Marquardt'); 结果还是一样。

于是使用1stOpt软件，网上只能下载到免安装的1.5版本，使用如下代码：
CODE:
Parameters Bt[10^10,10^12];
Parameters c[160,200];
Parameters m[-0.55,0.54];
Parameters n[-0.52,0.54];
Parameters p[0.65,1];
Parameters Rl;
Constant a=-150;
Constant b=150;
Constant Bx=-48030; //Bx,By,Bz是变参数，从文件中读出一行再解一次方程，对于每一组Bx,By,Bz需要求一次方程组
Constant By=29370;
Constant Bz=-19077.3;
//BatchFileModel;
//VarConstant Bx="Sheet1[A1:A240]";
//VarConstant By="Sheet1[B1:B240]";
//VarConstant Bz="Sheet1[C1:C240]";
//DataFile "D:\Mag1230.xls[Sheet1[A1:A240]]";
//RowDataSet;
     //Bx="Sheet1[A1:A240]";
     //By="Sheet1[B1:B240]";
     //Bz="Sheet1[C1:C240]";
//EndRowDataSet;
Function Bt*(3*(m*a+n*b+p*c)*a/Rl^5-m/Rl^3)/Bx-1=0;
Bt*(3*(m*a+n*b+p*c)*b/Rl^5-n/Rl^3)/By-1=0;
Bt*(3*(m*a+n*b+p*c)*c/Rl^5-p/Rl^3)/Bz-1=0;
m^2+n^2+p^2-1=0;
Rl-sqrt(a^2+b^2+c^2)=0;

结果为：
====== 结果 ======


迭代数: 148
计算用时(时:分:秒:毫秒): 00:00:00:482
计算中止原因: 达到收敛判定标准
优化算法: 最大继承法
函数表达式 1: bt*(3*(m*(-150)+n*150+p*c)*(-150)/rl^5-m/rl^3)/(-48030)-1-(0)
         2: bt*(3*(m*(-150)+n*150+p*c)*150/rl^5-n/rl^3)/29370-1-(0)
         3: bt*(3*(m*(-150)+n*150+p*c)*c/rl^5-p/rl^3)/(-19077.3)-1-(0)
         4: m^2+n^2+p^2-1-(0)
         5: rl-sqrt((-150)^2+150^2+c^2)-(0)
目标函数值: 0
Bt: 999812590445.556
c: 161.888089867363
m: 0.401079546362512
n: -0.0464426102549249
p: 0.914865171182263
Rl: 266.847809885904


====== 计算结束 ======
当然每次运行结果都会有点不同，但大致都在这个范围内。


因为我要继续一步作分析，Bx， By， Bz是测试数据，要从表格中读出的，对应每一组数据求一次方程组，相当于变常量，按照上面注释中采用的两种方法，都没有达到我预期的效果，软件把它也当做变量去求解了，也有可能是用的方法不对，但我从手册中找不到其它方法了。因为Matlab编程方便，可以实现我的需求，于是我又转回Matlab，用1stOpt解出非线性方程组的解作为初值用Matlab求解，依然求不出解，除了算法的设置，误差和迭代次数的设置都是因为解不出正确的解，自己调的，但是依然没用。
结果为：
====== 结果 ======


迭代数: 148
计算用时(时:分:秒:毫秒): 00:00:00:482
计算中止原因: 达到收敛判定标准
优化算法: 最大继承法
函数表达式 1: bt*(3*(m*(-150)+n*150+p*c)*(-150)/rl^5-m/rl^3)/(-48030)-1-(0)
         2: bt*(3*(m*(-150)+n*150+p*c)*150/rl^5-n/rl^3)/29370-1-(0)
         3: bt*(3*(m*(-150)+n*150+p*c)*c/rl^5-p/rl^3)/(-19077.3)-1-(0)
         4: m^2+n^2+p^2-1-(0)
         5: rl-sqrt((-150)^2+150^2+c^2)-(0)
目标函数值: 0
Bt: 999812590445.556
c: 161.888089867363
m: 0.401079546362512
n: -0.0464426102549249
p: 0.914865171182263
Rl: 266.847809885904


====== 计算结束 ======
当然每次运行结果都会有点不同，但大致都在这个范围内。


因为我要继续一步作分析，Bx， By， Bz是测试数据，要从表格中读出的，对应每一组数据求一次方程组，相当于变常量，按照上面注释中采用的两种方法，都没有达到我预期的效果，软件把它也当做变量去求解了，也有可能是用的方法不对，但我从手册中找不到其它方法了。因为Matlab编程方便，可以实现我的需求，于是我又转回Matlab，用1stOpt解出非线性方程组的解作为初值用Matlab求解，依然求不出解，除了算法的设置，误差和迭代次数的设置都是因为解不出正确的解，自己调的，但是依然没用。
CODE:
x0=[999999999999 161.9 0.401 -0.0465 0.9148 266.856];
options=optimset('Display','iter','MaxFunEvals',10^5,'TolFun',1e-18,'TolX',1e-18,'MaxIter',1e4,'Algorithm','Levenberg-Marquardt');
[x val]= fsolve(@solveMagneticDipoleEquation,x0,options);

请问：如何用Matlab可以解出这个方程组的解，或者用1stOpt可以实现我的分析需求？
有谁能给出一点建议，哪怕不一定是解决方法，Matlab里如何选择初值，用1stOpt解出的初值，为什么还是解不出来？
回复此楼

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : Mag1230.xls
    • 2015-03-02 16:14:15, 21.5 K

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-03-02 16:14:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mjjackey

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 月只蓝 at 2015-03-02 17:15:39
原方程组未知数之间的数量级相差过大，在不对原方程进行处理（诸如取对数降数量级之类）的情况下，求解的结果不稳定是正常的，而且此时方程用MALTAB求解比较困难，就是初值难以选取。建议用1stopt软件求解即可。
楼 ...

我也想完全用1stopt，关键是我不止要解出方程组的解，还要：“要继续一步作分析，Bx， By， Bz是测试数据，要从表格中读出的，对应每一组数据求一次方程组，相当于变常量，按照上面注释中采用的两种方法，都没有达到我预期的效果，软件把它也当做变量去求解了，也有可能是用的方法不对，但我从手册中找不到其它方法了。”
BatchFileModel;
VarConstant Bx="Sheet1[A1:A240]";
VarConstant By="Sheet1[B1:B240]";
VarConstant Bz="Sheet1[C1:C240]";

RowDataSet;
     Bx="Sheet1[A1:A240]";
     By="Sheet1[B1:B240]";
     Bz="Sheet1[C1:C240]";
EndRowDataSet;

我只有从网上下载的1.5破解版。还有我用最大继承法，目标函数值倒是为零，用你说的“标准简面体爬山法 + 通用全局优化法”，反倒收敛没它好。
赞一下 回复此楼
3楼2015-03-03 11:11:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 5 个回答

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与，应助指数 +1
mjjackey: 金币+5, ★★★很有帮助 2015-03-10 14:51:52
原方程组未知数之间的数量级相差过大，在不对原方程进行处理（诸如取对数降数量级之类）的情况下，求解的结果不稳定是正常的，而且此时方程用MALTAB求解比较困难，就是初值难以选取。建议用1stopt软件求解即可。
楼主编写的1stop代码略冗长，已修改如下：
CODE:
Constant a=-150, b=150, Bx=-48030, By=29370,Bz=-19077.3 ;
Parameters Bt[10^10,10^12],c[160,200],m[-0.55,0.54],n[-0.52,0.54],p[0.65,1]   ;

Function Bt*(3*(m*a+n*b+p*c)*a/(sqrt(a^2+b^2+c^2))^5-m/(sqrt(a^2+b^2+c^2))^3)/Bx-1=0;
Bt*(3*(m*a+n*b+p*c)*b/(sqrt(a^2+b^2+c^2))^5-n/(sqrt(a^2+b^2+c^2))^3)/By-1=0;
Bt*(3*(m*a+n*b+p*c)*c/(sqrt(a^2+b^2+c^2))^5-p/(sqrt(a^2+b^2+c^2))^3)/Bz-1=0;
m^2+n^2+p^2-1=0;

用简面体爬山法求解，效果良好：

化算法: 标准简面体爬山法 + 通用全局优化法
函数表达式 1: bt*(3*(m*(-150)+n*150+p*c)*(-150)/(sqrt((-150)^2+150^2+c^2))^5-m/(sqrt((-150)^2+150^2+c^2))^3)/(-48030)
            -1-(0)
         2: bt*(3*(m*(-150)+n*150+p*c)*150/(sqrt((-150)^2+150^2+c^2))^5-n/(sqrt((-150)^2+150^2+c^2))^3)/29370-1
            -(0)
         3: bt*(3*(m*(-150)+n*150+p*c)*c/(sqrt((-150)^2+150^2+c^2))^5-p/(sqrt((-150)^2+150^2+c^2))^3)/(-19077.3)
            -1-(0)
         4: m^2+n^2+p^2-1-(0)
目标函数值: 0
bt: 988812416631.896
c: 160.516748048499
m: 0.396983954369934
n: -0.0417362655864018
p: 0.916876122552939


MATLAB代码：
CODE:
function solve_nonl_eqs
clear all;clc
format long

x0=[1e12 160 0.4 -0.04  0.9];
options=optimset('Display','iter','Algorithm','Levenberg-Marquardt','TolFun',1e-21,'TolX',1e-21);
[x fval]=fsolve(@objfun,x0,options')


function  f=objfun(x)

%已知数
a=-150;
b=150;
Bx=-48030;
By=29370;
Bz=-19077.3;
%未知数
Bt=x(1);
c=x(2);
m=x(3);
n=x(4);
p=x(5);
Rl=sqrt(a.^2+b.^2+c.^2);



f(1)=Bt*(3*(m*a+n*b+p*c)*a/Rl.^5-m/Rl.^3)/Bx-1;
f(2)=Bt*(3*(m*a+n*b+p*c)*b/Rl.^5-n/Rl.^3)/By-1;
f(3)=Bt*(3*(m*a+n*b+p*c)*c/Rl.^5-p/Rl.^3)/Bz-1;
f(4)=m.^2+n.^2+p.^2-1;
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
2楼2015-03-02 17:15:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

引用回帖:
3楼: Originally posted by mjjackey at 2015-03-03 11:11:16
我也想完全用1stopt，关键是我不止要解出方程组的解，还要：“要继续一步作分析，Bx， By， Bz是测试数据，要从表格中读出的，对应每一组数据求一次方程组，相当于变常量，按照上面注释中采用的两种方法，都没有达 ...

1stopt软件及其教程：
http://muchong.com/bbs/viewthread.php?tid=4586018&fpage=1&target=blank
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
4楼2015-03-03 11:13:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mjjackey

新虫 (初入文坛)
引用回帖:
4楼: Originally posted by 月只蓝 at 2015-03-03 11:13:29
1stopt软件及其教程：
http://muchong.com/bbs/viewthread.php?tid=4586018&fpage=1&target=blank...

谢谢！这个软件的版本和我用的是一样的，手册我也下到了，网上广为流传的1.5版的破解版，太老了，一些功能都受限。
赞一下 回复此楼
5楼2015-03-03 12:14:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 5 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭