关于超越方程组的fsolve解法问题

小妮妮

我用fsolve函数求解超越方程组，可以求解出数值解，但是得到的结果和文献上的结果是有一定差别的，不知道问题出现在那里，请高手帮忙看看，万分感谢！下附我的程序。
%%%%%%%%%%%%下面是mainprogram%%%%%%%%%%%%%%%%%%%%%%
clear all;clc;close all;
%load ('m1.dat')
global js;
global q;
Tc =[];
q_par  = [0.9 1.28 2];
js_par = 2:-0.01:0.0;
n0     = [0.1 0.1 5.5];
for q = q_par
    js = js_par(1);
    n(1,:)=fsolve(@myfun1,n0,optimset('Display','off'));
    for i = 2:length(js_par)
        js = js_par(i);
        n(i,:)=fsolve(@myfun,n(i-1,:),optimset('Display','off'));
    end
    Tc = [Tc n(:,3)]
end
plot(Tc,js_par);
str =[];
for i=1:length(q_par)
    str = strvcat(str,num2str(q_par(i)));
end
legend(str,4);
xlabel('Tc/Ja')
ylabel('Js/Ja')
figure;
plot(Tc,js_par);
legend(str,4);
xlabel('Tc/Ja')
ylabel('Js/Ja')
%axis([0.08 2 0 2]);
%%%%%%%%%%%%以下是求解的超越方程组%%%%%%%%%%%%%%%%%%%%
function y=myfun(n)
global js
global q
syms t sx1 sx2
je=1.0;
qs=2.0;
k=1.0;
sx1=n(1);
sx2=n(2);

t=n(3);
w31=sqrt(qs^2+4*js*qs*sx1+je*q*sx2);
w32=sqrt(q^2+4*q*sx2+je*qs*sx1+je*qs*sx1);
y(1)=qs/(2*w31)*tanh(w31/(2*k*t))-sx1;
y(2)=q/(2*w32)*tanh(w32/(2*k*t))-sx2;

N=[2*w31*coth(w31/(2*k*t))-4*js, -je, 0;
    -je, 2*w32*coth(w32/(2*k*t))-4, -je;
    0, -je, 2*w31*coth(w31/(2*k*t))-4*js];
y(3)=det(N);

[ 本帖最后由 eight 于 2007-8-15 20:47 编辑 ]

小妮妮

怎么没有人帮忙呢，弄了很长时间，不知道问题到底在那里，请高手帮忙看看吧，太着急了！:'(

rocwoods

现在比较忙，没有时间细考虑。给你推荐两个帖子：
http://forum.vibunion.com/forum/thread-48384-1-5.html
http://www.simwe.com/forum/thread-787570-1-1.html

小妮妮

楼上发的两个网址我看过了，我觉得我的问题不是初值选择的问题，还是请高手帮忙看看吧，万分感谢，现在课题全卡住了。:'(

rocwoods

请把你的方程公式、条件用图片形式或者word文档发上来，这样看太麻烦了。

小妮妮

我把公式和图形上传了，请高手帮帮忙，万分感谢！

[ 本帖最后由 eight 于 2007-8-20 20:38 编辑 ]

小妮妮

怎么没有人帮忙呢，太着急了，请高手帮忙看看！！

xjzuo

看了一下你的代码和公式,发现有几个问题:
1. myfun1, myfun 不同----是否是笔误?
2. w32似乎和公式不同, 是否是公式没编辑好?
你自己再确认一下,看看是否是问题所在.

小妮妮

1。n(1,:)=fsolve(@myfun1,n0,optimset('Display','off'));红色位置的内容应该为myfun是笔误
2。sx1,sx2,tc是未知数，myfun 中的w32没有问题 ---------你的公式中是x1,没有问题？－－－by xjzuo
谢谢楼上，麻烦再帮忙找找问题,在线等！！急！！！

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

小妮妮

w32的表达式见下面的附件

rocwoods

你把'display'后的'off'改称'on'会发现循环过程中有很多
Maximum number of function evaluations reached:
increase options.MaxFunEvals.
提示，这表示没有收敛。
有些是Optimizer appears to be converging to a point which is not a root.
Relative function value changing by less than max(options.TolFun^2,eps) but
sum-of-squares of function values is greater than or equal to sqrt(options.TolFun)
这表示根本不会收敛到方程的根。
只有提示是Optimization terminated: first-order optimality is less than options.TolFun.才表示那次求解成功了。改改初值吧，直到所有的提示都是Optimization terminated: first-order optimality is less than options.TolFun.才行
还是建议你研究研究那些帖子。特别是初值和optimset的设置。
告诉你个选初值技巧，根据成功求解出结果的初值，总结各个初值之间的比例关系进行设置，这样求解成功的概率大些。自己再研究研究吧。

[ 本帖最后由 rocwoods 于 2007-8-23 11:19 编辑 ]

无水1324

确实是不错的方法，方程有时候有多个解，在满足收敛的条件下，要求另外一个解怎么改变初值，有没有比较好的建议？谢谢！

rocwoods

我一般都是在整个变量域内均匀随机设置初值然后努力寻找有效的迭代方法。譬如随机行走法就是比较有效的迭代方法。非线性方程组可以转化为无约束全局最优化
这里有我以前的一个讨论。http://www.simwe.com/forum/viewthread.php?tid=791513

xjzuo

再次看了一下,你的代码好象并没有什么问题;运行结果也基本合理,不知你说的差别是指什么?-----是左边的那条竖线?

小妮妮

1。当q_par=2.0时，文献上的曲线与横轴（图形的上边框）和纵轴的截距分别为1.5和1.33左右，而我画出来的图形与横轴和纵轴的截距分别为1.7和1.1左右，在Tc比较小的时候得到的解不理想，在程序中显示exitflag 发现小于零，说明没有收敛，怎么调试？
2。当q_par=1.065时，画出来的图形和文献上的图形状上就不相同，与图形的上边框和图形的下边框的截距分别为1.78和0.69左右，而文献上大约为0.27和1.56，用程序运行的结果都符合收敛条件。
3。当q_par=0.9时，画出来的图形和文献上的图形状上相似，与图形的上边框和图形的下边框的截距分别为1.78和0.79左右，而文献上大约为0.6和1.56，用程序运行的结果都符合收敛条件。
      此外，对于初值得选择问题，sx1和sx2得值都在0.5以下，t的值就是用从文献图中得到的值作为初值，具体如下
当q_par=2.0时，n0= [ unifrnd(0,0.5)  unifrnd(0,0.5)  1.5]
当q_par=1.065时，n0= [ unifrnd(0,0.5)  unifrnd(0,0.5)  1.56]
当q_par=0.9时，n0= [ unifrnd(0,0.5)  unifrnd(0,0.5)  1.56]
     但是始终都得不到理想的结果，本人实在不知问题出现在那里，也想了很长时间 ，一直解决不了了，现在课题全都卡住了，走投无路了，还是请高手帮忙看看，万分感激！！！！急急急！！！