求教：非线性方程组迭代解法的编程问题

camelry

方程组如附件，X1=VA，X2=6A

按照书上的例子，先编写bro：
function [y,n]=bro(x0,eps)
%布鲁顿迭代法
if nargin==2, eps=1.0e-6
elseif nargin<2,  error, return; end
A=eye(length(x0)); x1=x0-fbro(x0)/A; j=1;
while (norm(x1-x0)>=eps) & (j<=1000)
    x0=x1; x1=x0-fbro(x0)/A;
    p=x1-x0; q=fbro(x1)-fbro(x0);
    A=A+(q-p*A)'*p/norm(p); j=j+1;
end
y=x1; n=j;

%% 再编写fbro
function y=fbro(x)
d=5; ve=3; ro=0.4;
y(1)=(x(1)*normcdf((log(x(1)/d)+0.02+0.5*x(2)*x(2))/x(2))-exp(-0.02)*d*normcdf((log(x(1)/d)+0.02-0.5*x(2)*x(2))/x(2)))-ve;
y(2)=(x(1)*x(2)*normcdf((log(x(1)/d)+0.02+0.5*x(2)*x(2))/x(2)))/ve-ro;
y=[y(1) y(2)];


计算不出来结果，请大侠指点一二。谢谢


[ 本帖最后由 ChaChing 于 2010-3-1 08:26 编辑 ]

camelry

再命令窗口输入>>[y,n]=bro([2 0.5],eps)

算出结果。。。。。
y =
  1.0e+003 *
    0.4800   -9.1507

n =
        1001


可是X2应该是大于零的才对。。。。。。

[y,n]=bro([1 0.5],eps)

y =
    8.3036    0.0006

n =
     1001

参数不同，答案也不同。。。。。很多时候参数输入，就报错


上面的也不是答案，是我J值设定的太小了，1000


方程里含正态分布，可以用牛顿法解麽？？

[ 本帖最后由 xinyuxf 于 2007-5-8 09:28 编辑 ]

camelry

请高手指点，用fsolve可以解决麽？

如果这个问题解决了，再变d/ve/ro为变量，读取一个2000*3的矩阵，求得x1/x2的2000个答案。

我用循环语句控制，i=i+1，可以求解，问题是我不会把答案输入，每次都是拷贝屏幕上的答案到TXT里整理半天，好傻哦，请高手指点如何保存答案？

*************************
我就这么整理，删除其他文字，只保留答案，再贴进EXCEL里面。。。。
va =

  1.0e+007 *

  Columns 1 through 5

                   0   0.000015191744873   0.000001959011473   0.000004423352757   0.000007518352076

  Columns 6 through 10

   0.000004005238744   0.000004556613173   0.000015579341581   0.000006500117041   0.000002330693525

  Columns 11 through 15

   0.000003743739380   0.000022860635561   0.000093524485905   0.000003339738823   0.000002639275663

  Columns 16 through 20

   0.000010126396430   0.000008367559064   0.000012712029811   0.000004602640151   0.000002300267011

  Columns 21 through 25

camelry

用牛顿、布鲁顿迭代，都出不来结果

换成fsolve,JIU :victory:

就出来了，太好了

******************
>> x0=[1;1];
>> options=optimset('Display','iter');
>> [x,fv]=fsolve(@myfun,x0,options)
                                         Norm of      First-order   Trust-region
Iteration  Func-count     f(x)          step         optimality    radius
     0          3         8.83889                         0.721               1
     1          6         5.56837              1           1.36               1
     2          9         2.13659            2.5          0.363             2.5
     3         10         2.13659           6.25          0.363            6.25
     4         13         1.41694         1.5625          0.499            1.56
     5         14         1.41694        3.90625          0.499            3.91
     6         17        0.905391       0.976562          0.225           0.977
     7         18        0.905391        2.44141          0.225            2.44
     8         21        0.633925       0.610352          0.493            0.61
     9         24        0.398121        1.52588          0.949            1.53
    10         27     0.000471548       0.546743         0.0569            1.53
    11         30    1.98947e-009     0.00828717      8.48e-005            1.53
    12         33    7.81081e-020    3.3332e-005      6.61e-010            1.53
Optimization terminated: first-order optimality is less than options.TolFun.
x =
    7.9008
    0.1520

fv =
  1.0e-009 *
    0.1175
   -0.2536
**************************

下面再接着研究如何带入2000*3的变量矩阵，呵呵

还有导出答案

camelry

function F=myfun(x)
d=5; ve=3; ro=0.4;
F=[((x(1)*normcdf((log(x(1)/d)+0.02+0.5*x(2)*x(2))/x(2))-exp(-0.02)*d*normcdf((log(x(1)/d)+0.02-0.5*x(2)*x(2))/x(2)))-ve);((x(1)*x(2)*normcdf((log(x(1)/d)+0.02+0.5*x(2)*x(2))/x(2)))/ve-ro)];

[ 本帖最后由 ChaChing 于 2010-3-1 08:32 编辑 ]

camelry

搞定了第二步，可以带入变量了

>> mymain
Optimization terminated: first-order optimality is less than options.TolFun.
Optimization terminated: first-order optimality is less than options.TolFun.  (repeat)

ans =

   39.4965   13.7142   27.6819   30.1427  146.0034   15.6690   29.6886   43.6006   30.3374   54.3623   48.7404

>>

[ 本帖最后由 ChaChing 于 2010-3-1 08:27 编辑 ]

camelry

function [x1,x2,dd]=mymain
i=1;
x1=[];
x2=[];
dd=[];
ld=[32.878;10.915;23.233;28.777;141.119;12.671;24.729;39.577;28.41;52.58;46.987];
ro=[0.427;0.393;0.488;0.407;0.493;0.342;0.516;0.387;0.415;0.353;0.41];
ve=[7.274;3.016;4.919;1.937;7.708;3.249;5.466;4.809;2.492;2.824;2.686];
std=[32.878;10.915;18.96;28.777;136.364;12.671;23.027;38.354;21.049;52.482;46.725];
while i<=11
    x0=[1;1];
    x=fsolve(@myfun2,x0,'[]',ro(i),ve(i),ld(i));
    x1(i)=x(1);
    x2(i)=x(2);
    dd(i)=(x(1)-std(i))/(x(1)*x(2));
    i=i+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function F=myfun2(x,ro,ve,ld)
F=[((x(1)*normcdf((log(x(1)/ld)+0.02+0.5*x(2)*x(2))/x(2))-exp(-0.02)*ld*normcdf((log(x(1)/ld)+0.02-0.5*x(2)*x(2))/x(2)))-ve);((x(1)*x(2)*normcdf((log(x(1)/ld)+0.02+0.5*x(2)*x(2))/x(2)))/ve-ro)];

camelry

现在已经解决了参数的传递问题，最后一步，就是参数的导入，以及方程解的导出

这个该用什么命令呢？


为什么10楼的程序，会出来9楼的结果？？？只显示了X1，没有显示X2......

修改了options
options=optimset('Display','off');
    x=fsolve(@myfun2,x0,options,ro(i),ve(i),ld(i));

结果就没有前面的信息了，只有X1的解，可是X2呢？？

[ 本帖最后由 ChaChing 于 2010-3-1 13:36 编辑 ]

camelry

>> mymain

ans =

   39.4965   13.7142   27.6819   30.1427  146.0034   15.6690   29.6886   43.6006   30.3374   54.3623   48.7404


>> [x1,x2,dd]=mymain

%%%%%%%%%%%%%%%%%%%

hahahahah，自己全部搞定！～～～～～～～～～太佩服自己了，接触MATLAB20个小时，就解决了方案。

>> [x1,x2,dd]=mymain

x1 =

   39.4965   13.7142   27.6819   30.1427  146.0034   15.6690   29.6886   43.6006   30.3374   54.3623   48.7404


x2 =

    0.0790    0.0866    0.0877    0.0263    0.0266    0.0709    0.0965    0.0428    0.0343    0.0184    0.0227


dd =

    2.1214    2.3574    3.5911    1.7227    2.4851    2.6968    2.3246    2.8107    8.9279    1.8828    1.8188

:victory: :victory: :victory:

[ 本帖最后由 xinyuxf 于 2007-5-8 09:23 编辑 ]

hunter_009

自导自演!

camelry

5.8就要交导师, 5.1长假期间，请教了周围的朋友，没有人会MATLAB
昨天被逼急了，上网找，搜索到了这个论坛，20点发贴，也没人帮忙
只好自己翻这个论坛的贴，足足翻了60页，几乎是一年的全部帖子了，
发现自己一开始用布鲁顿是错误的，改用fsolve，就OK了

不过MATLAB的基本的命令都不会，就是看着别人的程序，找葫芦画瓢
一个人，就这么摸索着，在漫漫黑夜里。。真的是黑夜啊，我昨夜通宵

终于自己解决了，难道不高兴么？ 楼上的说什么自编自演。呵呵。。

10楼的真狭隘。不过我一点也不生气。论文完成了，真是爽啊！

[ 本帖最后由 ChaChing 于 2010-3-1 08:31 编辑 ]

camelry

评分分数: 振动币 +14 个
操作理由: 虽然你的自说自话让人觉得...
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
我自説自话，是因为半夜三更的，没人来帮助我，我只好自己记录下来摸索的痕迹，万一不成功，白天有高手看见了可以指点一下。

例如，我现在还不知道如何导入大量的数据，只能事先手工把变量先输入M文件里。

camelry

例如，解出来的X1/X2/DD，我预想的是竖列的，结果答案里是横排的，拷贝回EXCEL出错，超过255列了，我不知道转置的命令。只好一批批的200个200个拷贝进EXCEL，呵呵，自己汗一个。。。。

xinyuxf

你摸索问题的精神很可嘉，如果你有问题或者找到了答案想要分享，都可以发在论坛上，这也正是论坛所提倡的。
    但是，我认为在论坛上发帖、回帖，内容应该都是你所要表达的思想的关键、要点，应该是比较简要、易懂的，而不是把论坛当作你解决问题这一过程的详细的记录本。你也提到了，自己看了60页的网页，肯定很累很繁琐，你可以换个角度思考，试想如果大家都像你这样发帖子，像是个演草本，把所有东西都罗列上去，其他人浏览帖子也会不舒服的。
    我想其他人说的话，没有什么恶意，只是可能表达方式不对。尤其是xjzuo，都已经给你加分，是对你努力的肯定，对你发帖存在的问题提出意见也是应该的，他是本版的版主，有责任指正在本版发帖会员的不当。

dingd

用1stOpt求解：

代码：
Constant d=5,ve=3,ro=0.4;
Function
(x1*normcdf((log(x1/d)+0.02+0.5*x2*x2)/x2)-exp(-0.02)*d*normcdf((log(x1/d)+0.02-0.5*x2*x2)/x2))-ve;
(x1*x2*normcdf((log(x1/d)+0.02+0.5*x2*x2)/x2))/ve-ro;

结果：
x1: 8.00023237371198
x2: 0.161444269349853