Matlab教学实验常见问题反馈(by dsj6700417)

建不了的模。

我讲一下大家常犯的错误：

1.有同学要拟合的目标函数是a*x+b，使用了最小二乘拟合lsqcurvefit，或者非线拟合nlinfit函数，
  应该使用多项式拟合polyfit.具体到此例使用polyfit(x,y,1)

2.有同学要数值求解多参数方程组，使用了fzero.
  fzero只能求单参数方程的零点.多参数方程的求解用fsolve
    试看下面的例子：
x=[0.0350;0.1793;0.5573];%初值
    options=optimset('MaxFunEvals',30000,'MaxIter',30000);
    [x,Fval,exitflag] = fzero(@fx,x,options);x
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
    function f=fx(x)
    F=x(1);K=x(2);K_a=x(3);
    f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
     -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
        -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

上面的程序格式上是正确的，但报错的原因在于错误使用了fzero，改为fsolve即可运行

3.我以lsqcurvefit为例，讲一下这个函数使用上的技巧，并解释大家的疑问
先讲使用lsqcurvefit的四个注意事项：
(1)注意模型函数的正确性，
(2)注意初值的选择
(3)注意optimset的设置
(4)注意res的大小，和exitflag的返回值

  我先用下面的代码分析一下：
    clear
    x=0:9;
    y=[1        2.31969        4.50856        6.90568        6.00512        5.56495        5.32807        7.56101        8.9392        9.5817];

    f=inline('1/(a(1)+(a(2)-a(1))*exp(-a(3)*x))','a','x');%logstic函数
    [xx ,res]=lsqcurvefit(f,[1 1 1],x,y)
复制代码
这是lsqcurvefit最常用的使用方法.上面程序的报错原因在于，inline函数里面的分母是个数组，数组除法要用点除，
在1后面加个点，即1./(a(1)+(a(2)-a(1))*exp(-a(3)*x))，程序就可以运行了.

修改后的程序，残差res很大，任凭你怎么设置optimset的精度及容差，或者修改初值，这个残差res都很大，这时需要考虑函数模型的正确性。logstic函数是个单调函数，而数据的散点图不是单调的，中间一段有下降部分，这个是导致res过大的原因，所以兔子问题我们根据logstic函数的单调性分为三段拟合。

我把问题提的尖锐一些：你不可能用指数函数exp(k^2*x),去拟合散点图呈正弦分布的数据，而要求残差很小；你不可能用一个恒正的函数模型去拟合全负的数据;你又不能用一个单调增函数去拟合一组单调减的数据。

后面这种情况，兔子问题的第二段曲线拟合就是个例子，初值都取正数时，logstic函数就是单调增函数，这样去拟合一组单调减的数据，结果是很差的，解决办法是改最后一个初值为负.如果你细心分析的话，会发现第三个初值为负就可以保证logstic函数是单调减的。

这个例子说明：
(1)要保证函数模型正确，不能使用错误的或者不合适的函数模型；
(2)构建函数时，要注意数组运算符的正确使用；
(3)初值的设置与函数模型的性质有关，要学会大体分析函数模型的性质.拟合残差过大时，试着多改几次初值；
(4)注意观察残差res的大小，res越小，拟合效果越好；
(如果使用[x,resnorm,residual,exitflag]=lsqcurvefit(f,[1 1 1],x,y) ，则exitflag返回-1时，说明拟合很成功)


我们看一下初值选取要注意的问题，先看下例：
    clear
    t=[ 3   18   20];
    y=[ 763.9000  76.3900  53.4000];
    f=inline('a(1)*0.1*(exp(-a(2)*t)-exp(-a(3)*t))/(a(3)-a(2))','a','t');

    ff=optimset;
    ff.MaxFunEvals=100000;
    ff.MaxIter=10000;
    [xx ,res]=lsqcurvefit(f,[1 ,1 ,1],t,y,[],[],ff)
复制代码

程序报错的原因在于inline函数的最后：(a(3)-a(2))在取初值时等于了0，导致分母为零，优化直接中止.一般情况下，大家习惯取初值为全1矩阵.这也要视具体函数而定.

我们把初值改为[.1,.1,.2]，程序即可运行，但此后不管我们怎么修改初值(要保证后两个不相等)，残差res都大的惊人，这时的原因就在于函数模型给的不对(不是指函数格式不正确，是指问题的函数模型有误)，或者说是函数模型使用的数据不对，上例正是这个原因：函数模型错误


再看下面的例子(注意 fsolve也属于优化函数)
x=[1 1 2];%初值
[x,Fval,exitflag] = fsolve(@fx,x);x
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
    function f=fx(x)
    F=x(1);K=x(2);K_a=x(3);
    f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
     -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
        -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码


上面初值的后面两个不相等也是为了避免函数分母为零.

为了得到优化结果,好多人都会用下面的办法：

随意设置一个初值，把运行结果作为初值，反复执行，直到初值和结果相等，即迭代至不动点为止。这个办法是很耗功耗时的，有一个一劳永逸的办法是：在optimset里面调大函数调用次数MaxFunEvals和迭代次数MaxIter，但很多情况上面的办法都不起效果，比如上面的程序改为：

  x=[1 1 2];%初值
    options=optimset('MaxFunEvals',30000,'MaxIter',30000);
    [x,Fval,exitflag] = fsolve(@fx,x,options);x
复制代码


不管你把MaxFunEvals和MaxIter的值调的再大，结果都不会收敛。这时必须更换初值的类型，比如一个好的运气，把初值改为

x=[3;4;200];%初值
    options=optimset('MaxFunEvals',30000,'MaxIter',30000);
    [x,Fval,exitflag] = fsolve(@fx,x,options);x
复制代码
优化就中止了.


总结一下，优化时综合考虑以下问题：
(1)函数模型的正确性，函数的正确书写
(2)初值的选择，主要与模型函数的性质有关
(3)optimset的设置，主要用于设置精度，容差，最大调用次数，最大迭代次数
(4)学会查看优化效果：res和exitflag

优化问题本来就是很复杂的，要想得到需要的结果需要设置很多选项和考虑很多因素，鉴于此，不提倡大家使用cftool!

我把血药浓度问题中，三个参数的求取程序附在后面，提供给大家一个参考，大家有兴趣的可以试着把自己有错误的程序，按照上面的错误警醒，一步步修改一下，体会一下优化需要注意的细节.
%%确定F,K,K_a的值
    x=[0.0350;0.1793;0.5573];%x=[3;4;200];%初值
    options=optimset('MaxFunEvals',30000,'MaxIter',30000);
    [x,Fval,exitflag] = fsolve(@fx,x,options);x   % 返回F,K,K_a的具体数值，x=[F,K,K_a]
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
        function f=fx(x)
        F=x(1);K=x(2);K_a=x(3);
        f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
         -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
            -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

深V礼

顶顶，好样的