求助，求解一个非线性方程组【已经编辑初值问题】

manking

我是matlab的菜鸟，没有用matlab算出来，希望哪位大侠帮忙求解一下
方程组中除x，y，m，n外均为已知量，R=50，D=10，那个像B的bai(这里是读音，打不出来这个字母，呵呵)的值为0.75，p0=10e6
方程组如图所示：

非常感谢各位兄弟参与讨论，我这几天准备整理一下，把结果及程序贴上来
关于matlab解，很多兄弟说初值的问题，我这里解释一下，m和n的值都是在大于0小于50的范围内，x，y和p0的值应该不会相差一个数量级，若解得的值不在上面范围之内，估计就是有问题了:loveliness:

我这里根据第一位兄弟提供的方法试验了一下结果如下：
"Title "feixianxingfangchengzu";
//Parameters x,y,m,n;
//Variable ;
Constant R=50, D=10, bta=0.75, p=10e6;
Function p+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*y/(R+D-n)^3=R^3*y/(R-n)^3;
         p+2*bta*R^3*x/(5*R+2*D-m)^3+2*bta*R^3*x/(R+D-m)^3+2*bta*R^3*y/(3*R+D-m)^3=R^3*y/(R+n)^3;
         p+2*bta*R^3*x/(5*R+3*D-m)^3+2*bta*R^3*y/(3*R+2*D-n)^3+2*bta*R^3*y/(R+D+n)^3=R^3*y/(R-m)^3;
         p+2*bta*R^3*x/(7*R+3*D-m)^3+2*bta*R^3*y/(5*R+2*D-n)^3+2*bta*R^3*y/(3*R+D+n)^3=R^3*x/(R+m)^3;

====== 结果 ======

迭代数: 32
计算用时(时:分:秒:毫秒): 00:00:07:63
计算中止原因: 达到收敛判定标准
优化算法: 包维尔法 + 通用全局优化法
函数表达式 1: 10000000+2*0.75*50^3*x/(3*50+2*10-m)^3+2*0.75*50^3*x/(3*50+2*10-m)^3+2*0.75*50^3*y/(50+10-n)^3-(50^3
            *y/(50-n)^3)
         2: 10000000+2*0.75*50^3*x/(5*50+2*10-m)^3+2*0.75*50^3*x/(50+10-m)^3+2*0.75*50^3*y/(3*50+10-m)^3-(50^3*y
            /(50+n)^3)
         3: 10000000+2*0.75*50^3*x/(5*50+3*10-m)^3+2*0.75*50^3*y/(3*50+2*10-n)^3+2*0.75*50^3*y/(50+10+n)^3-(50^3
            *y/(50-m)^3)
         4: 10000000+2*0.75*50^3*x/(7*50+3*10-m)^3+2*0.75*50^3*y/(5*50+2*10-n)^3+2*0.75*50^3*y/(3*50+10+n)^3-(50^3
            *x/(50+m)^3)
目标函数值: 3.72529029846191E-9
x: 11230950.1093155
m: -0.610358781362475
y: 32159900.0699049
n: 7.61505721423252

====== 计算结束 ======
m的值不在范围之内，等我用各位高手提供的matlab方法试验一下
看看 效果如何:loveliness:

[ 本帖最后由 xjzuo 于 2007-1-8 10:24 编辑 ]

dingd

1stOpt求解非线性方程组比Matlab来的方便、几乎不用编程，关键是不需猜初值（一件很头疼的事），下面是1stOpt代码和结果：

Constant R=50, D=10, bta=0.75, p=10e6;
Function p+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*y/(R+D-n)^3=R^3*y/(R-n)^3;
         p+2*bta*R^3*x/(5*R+2*D-m)^3+2*bta*R^3*x/(R+D-m)^3+2*bta*R^3*y/(3*R+D-m)^3=R^3*y/(R-n)^3;
         p+2*bta*R^3*x/(5*R+2*D-m)^3+2*bta*R^3*y/(3*R+2*D-n)^3+2*bta*R^3*y/(R+D+n)^3=R^3*y/(R-m)^3;
         p+2*bta*R^3*x/(7*R+3*D-m)^3+2*bta*R^3*y/(5*R+2*D-n)^3+2*bta*R^3*y/(3*R+D+n)^3=R^3*y/(R+m)^3;


x: 28352109.1808853
m: 3.70370471729419
y: 13271271.0409974
n: 15.5934304924452

manking

lxq, yangzj, xueyi, xjzuo, huright
几位版主能帮下么:@)

flybaly

function f=try111(t)
R=50;D=10;bta=0.75;p=10e6;
x=t(1);
y=t(2);
m=t(3);
n=t(4);
f1=p+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*y/(R+D-n)^3-R^3*y/(R-n)^3;
f2=p+2*bta*R^3*x/(5*R+2*D-m)^3+2*bta*R^3*x/(R+D-m)^3+2*bta*R^3*y/(3*R+D-m)^3-R^3*y/(R+n)^3;
f3=p+2*bta*R^3*x/(5*R+2*D-m)^3+2*bta*R^3*y/(3*R+2*D-n)^3+2*bta*R^3*y/(R+D+n)^3-R^3*y/(R-m)^3;
f4=p+2*bta*R^3*x/(7*R+3*D-m)^3+2*bta*R^3*y/(5*R+2*D-n)^3+2*bta*R^3*y/(3*R+D+n)^3-R^3*y/(R+m)^3;
f=[f1;f2;f3;f4];
在命令行中输入：
x0 = [-400.21545164; 13226445.6938747;3.818954534493;13.775603875429];% Make a starting guess at the solution
options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@try111,x0,options)  % Call optimizer
结果：

                                         Norm of      First-order   Trust-region
Iteration  Func-count     f(x)          step         optimality    radius
     0          5    1.84006e+013                     8.63e+011               1
     1         10     1.7908e+013              1      1.47e+011               1
     2         15     1.7874e+013            2.5      2.85e+010             2.5
     3         20    1.78732e+013           6.25       5.2e+009            6.25
     4         25    1.78731e+013         15.625      1.26e+009            15.6
     5         30    1.78727e+013        39.0625      2.44e+008            39.1
     6         35    1.78718e+013        97.6562      5.01e+007            97.7
     7         40    1.78696e+013        244.141      6.37e+007             244
     8         45    1.78642e+013        610.352      1.43e+008             610
     9         50    1.78506e+013        1525.88       3.6e+008       1.53e+003
    10         55    1.78165e+013         3814.7         9e+008       3.81e+003
    11         60    1.77316e+013        9536.74      2.25e+009       9.54e+003
    12         65      1.752e+013        23841.9      5.58e+009       2.38e+004
    13         70    1.69966e+013        59604.6      1.37e+010       5.96e+004
    14         75    1.57224e+013         149012      3.31e+010       1.49e+005
    15         80    1.27519e+013         372529      7.45e+010       3.73e+005
    16         85     6.6813e+012         931323      1.35e+011       9.31e+005
    17         90    1.78718e+010   2.32831e+006      1.91e+010       2.33e+006
    18         95         2910.14         126547       2.1e+007       5.82e+006
    19        100    1.07262e-010         49.771           10.7       5.82e+006
    20        105    3.46945e-018   5.83906e-006         0.0011       5.82e+006
    21        106    3.46945e-018    1.8823e-009         0.0011       5.82e+006
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).
x =
  1.0e+007 *
  -0.40075001228953
   1.32264457905741
   0.00000038189545
   0.00000137756039
fval =
  1.0e-008 *
                  0
                  0
                  0
   0.18626451492310
和二楼的结果不一样是正常的。因为二楼的公式（第二个公式的右侧分母）好像错了。
1stOPT中的结果为：
====== 结果 ======

迭代数: 56
计算用时(时:分:秒:毫秒): 00:00:02:344
计算中止原因: 达到收敛判定标准
优化算法: 麦夸特法(Levenberg-Marquardt) + 通用全局优化法
函数表达式 1: 10000000+2*0.75*50^3*x/(3*50+2*10-m)^3+2*0.75*50^3*x/(3*50+2*10-m)^3+2*0.75*50^3*y/(50+10-n)^3-50^3
            *y/(50-n)^3-(0)
         2: 10000000+2*0.75*50^3*x/(5*50+2*10-m)^3+2*0.75*50^3*x/(50+10-m)^3+2*0.75*50^3*y/(3*50+10-m)^3-50^3*y
            /(50+n)^3-(0)
         3: 10000000+2*0.75*50^3*x/(5*50+2*10-m)^3+2*0.75*50^3*y/(3*50+2*10-n)^3+2*0.75*50^3*y/(50+10+n)^3-50^3
            *y/(50-m)^3-(0)
         4: 10000000+2*0.75*50^3*x/(7*50+3*10-m)^3+2*0.75*50^3*y/(5*50+2*10-n)^3+2*0.75*50^3*y/(3*50+10+n)^3-50^3
            *y/(50+m)^3-(0)
目标函数值: 0.000331555493175983
x: -4007500.12264468
m: 3.818954534493
y: 13226445.7907672
n: 13.775603875429

====== 计算结束 ======

[ 本帖最后由 flybaly 于 2007-1-5 19:24 编辑 ]

xjzuo

对于这种"买卖贴"本想直接删除的...
另: 这种问题自己动动手就可解决, 不鼓励帮其写程序计算.

flybaly

原帖由 xjzuo 于 2007-1-5 19:42 发表
对于这种"买卖贴"本想直接删除的...
另: 这种问题自己动动手就可解决, 不鼓励帮其写程序计算.

本意是指正二楼的错误之处。
既然程序写出来了，也就顺便贴上，不只是为了帮楼主解答，也是借花献佛。
正好比较一下MATLAB和1STOPT解的差别，也为看本贴的人不白看。
来本论坛都是为了学习，不必为新来人这么苛刻。
个人认为看程序学习还是比较有效的一种。


%%%---------------------------------------------------%%%
你愿意帮他这没有什么问题,我也看出了你想比较一下两种软件的计算结果.
我只是想强调一点,对于此类买卖贴,不鼓励回帖.
对于一般正常讨论贴,我本身也回过.论坛上也可搜到fsolve贴．
我本人也没有对新来人苛刻的意思．
%%%------------------------------------------------------------%%%
By xjzuo

[ 本帖最后由 xjzuo 于 2007-1-5 23:27 编辑 ]

manking

原帖由 xjzuo 于 2007-1-5 19:42 发表
对于这种"买卖贴"本想直接删除的...
另: 这种问题自己动动手就可解决, 不鼓励帮其写程序计算.

:@L 这位兄弟谅解
本来自己准备动手做的
因为有很多事情就耽搁来的
没有买卖的意思:@(
只是为了表示谢意阿:handshake

manking

原帖由 flybaly 于 2007-1-5 19:19 发表
function f=try111(t)
R=50;D=10;bta=0.75;p=10e6;
x=t(1);
y=t(2);
m=t(3);
n=t(4);
f1=p+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*x/(3*R+2*D-m)^3+2*bta*R^3*y/(R+D-n)^3-R^3*y/(R-n)^3;
f2=p+2*bta*R^3 ...

谢谢兄弟指教:handshake
请兄弟任意选一对
我会pm你密码的:loveliness:

509562379
509562380

509562575
509562576

509562177
509562178

manking

原帖由 flybaly 于 2007-1-5 20:08 发表


本意是指正二楼的错误之处。
既然程序写出来了，也就顺便贴上，不只是为了帮楼主解答，也是借花献佛。
正好比较一下MATLAB和1STOPT解的差别，也为看本贴的人不白看。
来本论坛都是为了学习，不必为新来人 ...

:loveliness: 再次感谢兄弟:handshake

xjzuo

除了上面对flybaly的答复,我想强调一点: 对于一般简单的问题,
直接写出程序不会比"先给思路,让其先动手,然后再讨论"效果好.
对于复杂的问题,如果有时间参与讨论,当然鼓励直接写出自己的见解.
同时也希望flybaly不要误解我的本意.

[ 本帖最后由 xjzuo 于 2007-1-5 23:39 编辑 ]

flybaly

原帖由 xjzuo 于 2007-1-5 23:37 发表
除了上面对flybaly的答复,我想强调一点: 对于一般简单的问题,
直接写出程序不会比"先给思路,让其先动手,然后再讨论"效果好.
对于复杂的问题,如果有时间参与讨论,当然鼓励直接写出自己的见解.
同时 ...

还是讨论问题吧。
这本没有谁对谁错的问题，只要目的正确就好。个人还是感觉这方面还是仿真论坛做得好一些。

manking

引用xjzuo

对于这种"买卖贴"本想直接删除的...
另: 这种问题自己动动手就可解决, 不鼓励帮其写程序计算.

也许由于是新人，不是很熟悉本论坛的的一些习俗，不对之处还望见谅
以后不会出现兄弟所说的这种发贴方法及方式
大大小小的论坛去的多了
没有规矩不成方圆，我已经编辑求助内容了
希望大家可以更好的继续讨论
当然偶说话还是要兑现的
有想需要号的朋兄弟直接pm我吧，我会第一时间pm你密码

引用flybaly 兄弟的话

本意是指正二楼的错误之处。
既然程序写出来了，也就顺便贴上，不只是为了帮楼主解答，也是借花献佛。
正好比较一下MATLAB和1STOPT解的差别，也为看本贴的人不白看。
来本论坛都是为了学习，不必为新来人这么苛刻。
个人认为看程序学习还是比较有效的一种。

谢谢兄弟的理解与支持
我想也正是有兄弟这样的版友，论坛的交流与学习的气氛才会更好，更融通，论坛也会发展更好

这里也谢谢dingd 兄弟的的一块玉，希望大家本着学习的目的共同继续讨论、探讨这个问题，也给新手们一个全面了解的机会
:loveliness:

dingd

很抱歉啊，公式里弄错一个符号导致错误结果。

这类非线性方程组，Matlab可以解，但难点是初值不好定。五楼flybaly的Matlab代码中的初值不知是如何定出的？有何诀窍？不会是先参照 1stOpt的计算结果吧？:lol 一般人很难一下给出如此接近的初值。

用1stOpt（2.0更强），换一下算法，精度会更高一些。

xjzuo

原帖由 flybaly 于 2007-1-6 16:14 发表


还是讨论问题吧。
这本没有谁对谁错的问题，只要目的正确就好。个人还是感觉这方面还是仿真论坛做得好一些。

我之所以耐心解释,是因为我很想创造一个良好的讨论气氛,一个良好的讨论秩序,
要记住这里是非赢利性的学术讨论版,不是交易版,发贴还是请先注意一下.

另外,如果有些人不但不听管理, 还以为自己有理了,而且还作一些空洞的评论,那我想这里可能不适合你．
当然，我们会不断学习别的论坛，尤其是bainhome等领导的仿真论坛Matlab版，正是因为他们的管理有方，
对于垃圾贴的适当处理，才使得不少讨论贴都很有水准．

我之所以多次解释，也是为了进一步管理好论坛，希望论坛得到很好的发展，
使论坛不至于沦为乱七八糟的地方，使论坛讨论贴的水平不断得到提高．

同时我也相信，本论坛也有很多高手，例如本版块的Happy, eight等，
有了他们，我相信本论坛会不断发展，真正成为一个良好的学术讨论地．

最后再强调一遍：对于简单问题，建议先给思路，或先搜索论坛，还不会再来讨论；
对于复杂问题，有时间参与讨论的版友，鼓励直接发表自己的意见．尤其欢迎原创见解．

[ 本帖最后由 xjzuo 于 2007-1-6 22:13 编辑 ]

dingd

这道题有点意思啊，用1stOpt又找到一组解，不知对否，大家验证一下：

x: -6284980.20782905
m: 12.0166413640221
y: 21836463.6740884
n: 122.909601202379