matlab中fsolve默认终止误差,Matlab用fsolve函数求解非线性方程组不收敛，初值如何选择？...

%MATLAB编程求拉格朗日乘数法极大值过程如下

clear

syms L hc sigmaw sigmac w w0 h a a0 wt x ftmax b gf P CMOD G delta  M M1 Eq1 Phi Lamda l1 l2 l3 s W   %定义符号变量

%h=200;a0=60;b=100;l=800;gf=0.11;ftmax=2.55;W=392;%给各个变量赋值

w0=2*gf/ftmax;                          %w0为粘聚应力为零时的裂缝张开位移

sigmaw=ftmax*(w0-w)/w0;                 %使用单线性粘聚力模型求得的粘聚应力

w=(x-hc)*wt/(a-a0);                     %断裂过程区虚拟裂缝张开位移w沿界面高度的分布

G=int(int(sigmaw,0,w)*b,hc,hc+a-a0);    %双重积分求得整个截面的断裂能G=(b*ftmax*wt*(6*gf - ftmax*wt)*(a - a0))/(12*gf)

sigmac=ftmax*(h-a-hc)/hc;               %根据平截面假定，求得梁上边缘处混凝土的压应力sigmac与裂缝尖端极限抗拉强度ftmax的关系

sigmaw=ftmax*(w0-w)/w0;

Eq1=0.5*sigmac*b*(h-a-hc)-0.5*ftmax*b*hc-int(sigmaw*b,hc,hc+a-a0);

%hc=solve(Eq1,hc);                        %求解hc

hc =-(w0*a^2 - 2*w0*a*h + w0*h^2)/(2*a0*w0 + a*wt - a0*wt - 2*h*w0);

%求得hc的表达式，用a，wt，gf来表示。相当于消去了hc

%根据几何关系求得delta和CMOD以及wt的关系

% 0.5*P*delta=P*l*wt/8*(a-a0);

%M1=- (L*wt*(W/2 - (4*M)/L))/(8*a - 8*a0) - (b*ftmax*wt*(3*w0 - wt)*(a -a0))/(6*w0)

%建立了a和wt的平衡关系式M1=- (L*wt*(W/2 - (4*M)/L))/(8*a - 8*a0) - (b*ftmax*wt*(3*w0-wt)*(a - a0))/(6*w0)

hc =-(w0*a^2 - 2*w0*a*h + w0*h^2)/(2*a0*w0 + a*wt - a0*wt - 2*h*w0);

w=(x-hc)*wt/(a-a0);

sigmaw=ftmax*(w0-w)/w0;

M=ftmax*b*(h-a-hc)^3/(3*hc)+(1/3)*ftmax*b*hc^2+int(sigmaw*b*x,hc,hc+a-a0);%开裂截面弯矩M表达式,很长，暂不列出

P=4*M/L-0.5*W;           %p和M的关系式

M1=P*L*wt/(8*a-8*a0)-G;

%应用拉格朗日乘数法，可以建立一个拉格朗日函数Phi，将平衡方程M1左端乘以待求常数Lamda，并与M的表达式相加

%用拉格朗日乘数法可以得到一个关于a、wt、Lamda的三元非线性方程组

%用MATLAB求解可以得到临界有效裂缝长度ac、初始裂缝尖端的张口位移wt和开裂界面上最大弯矩Mmax

%最大断裂荷载可以通过下式求得   Pmax=4*Mmax/l-0.5*W       W是混凝土三点弯曲梁的自重

%可以看到，极限荷载Pmax主要取决于裂缝尖端极限抗拉强度ftmax和裂缝稳定扩展区域局部断裂能gf

Phi=M+Lamda*M1;                               %建立一个拉格朗日函数Phi

l1 =diff(Phi,a);                              %拉格朗日函数Phi对a求偏导数

l2 =diff(Phi,wt);                            %拉格朗日函数Phi对wt求偏导数

l3 =diff(Phi,Lamda);                          %拉格朗日函数Phi对Lamda求偏导数

%  h=200;a0=60;b=100;L=800;gf=0.11;ftmax=2.55;W=392;           % 参考冷科治论文数据

% h=203;a0=77.749;b=76;L=762;gf=0.14;ftmax=4.52;W=288;            % B系列梁B1

%  h=305;a0=60;b=76;L=1143;gf=0.14;ftmax=4.57;W=649;            % C系列梁

%给t0赋初值，再调用fsolve函数求解非线性方程组，如何收敛？？？Lamda的初值怎样选？？？？？？

t0=[92,12.5,50]

%t = fsolve('answer',t0)

options=optimset('Display','iter');

[t,fval] = fsolve(@answer,t0,options) %fval 是求解误差

输入上面几行后得到的运行结果如下，求解器提前结束了，并没有得到收敛解。郁闷中，如何调整初值，希望诸位大神不吝赐教。

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,

options.MaxFunEvals = 300 (the default value).

t =

102.4856   18.0198  333.4458

fval =

1.0e+003 *

4.9629

-0.2785

0.3396

下面是编写的函数

function f=answer(t)

%  h=200;a0=60;b=100;L=800;gf=0.11;ftmax=2.55;W=392;           % 参考冷科治论文数据

h=203;a0=77.749;b=76;L=762;gf=0.14;ftmax=4.52;W=288;            % B系列梁

%  h=305;a0=60;b=76;L=1143;gf=0.14;ftmax=4.57;W=649;            % C系列梁

a=t(1);wt=t(2)/1000;Lamda=t(3);

l1 =Lamda*((wt*((4*b*ftmax*wt*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)^2)/(3*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^3) - (4*b*ftmax*wt*(h - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^3)/(3*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)) - 4*b*ftmax*(a0 - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))*((wt*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax))/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^2 - ((4*a*gf)/ftmax - (4*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax) + 1) - (4*b*ftmax*((4*a*gf)/ftmax - (4*gf*h)/ftmax)*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax))/(3*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^2) - (2*b*ftmax^2*wt*(a0 - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^3)/(3*gf*(a - a0)^2) + (4*b*ftmax*((4*a*gf)/ftmax - (4*gf*h)/ftmax)*(h - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^3*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))/(3*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)^2) + (4*b*ftmax*(h - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^2*((wt*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax))/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^2 - ((4*a*gf)/ftmax - (4*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax) + 1)*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))/((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax) + (2*b*ftmax*h^2*wt*(a0 - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^2)/((a - a0)^2*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)) - (2*a^2*b*ftmax*wt*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)^2)/((a - a0)^2*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^3) - (2*b*ftmax*h^2*wt*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)^2)/((a - a0)^2*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^3) + (2*a^2*b*ftmax*wt^2*(a0 - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^2)/((a - a0)*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^2) + (2*b*ftmax*h^2*wt^2*(a0 - a + ((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)/(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax))^2)/((a - a0)*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^2) - (6*a^2*b*ftmax*wt^2*((2*a^2*gf)/ftmax + (2*gf*h^2)/ftmax - (4*a*gf*h)/ftmax)^2)/((a - a0)*(a*wt - a0*wt + (4*a0*gf)/ftmax - (4*gf*h)/ftmax)^4) + (2*b*ftmax^