Matlab--用最小二乘法确定指数函数的系数(涉及牛顿法解非线性方程组）

题目

最小二乘法解题步骤

显然，该方程组为非线性方程组，考虑用牛顿法求解。

牛顿法解非线性方程组

选取牛顿迭代的初值

牛顿法只具有局部收敛性，因此初始值的选取很重要。观察题中给出的数据，发现原函数 y=ae^{bx},有递增趋势,y值都较小且大于0，可以判断a,b的值应该同时为正且小于1，所以在区间(0,1)内随机选取初始值。又因为随机选取的初始值不能保证结果一定收敛，所以随机选取多组初始值进行多次计算，选取其中均方误差最小的一组结果，也即所有结果中效果最好的一组a,b作为最终结果。

计算结果

在区间(0,1)内随机选取50组初始值，其中能收敛且均方误差最小的一组是：
将所得函数和原始数据进行拟合比较（红线是对数据表中数据描点连线得到）

matlab代码

初始值是随机生成的，所以可能每次运行代码得到的结果都不一样

clear

%最小二乘法

x=1:19;

y=[0.898 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 4.76 6.46 6.53 10.9 16.5 22.5
35.7 50.6 61.6 81.8];

%建立目标函数

syms a b

F=0;

for i=1:length(x)

    F=F+(a*exp(b*x(i))-y(i))^2;

end

%偏导

Fa=diff(F,a);

Fb=diff(F,b);

%令偏导等于0，用牛顿法求解非线性方程组Fab

Fab=[Fa,Fb];

dFab(1,1)=diff(Fa,a);

dFab(1,2)=diff(Fa,b);

dFab(2,1)=diff(Fb,a);

dFab(2,2)=diff(Fb,b);

var=[a b];

delta=0.01;

t=0;

K(1,50)=0;%预分配内存

V(50,2)=0;

for i=1:50

v0=rand(1,2);%随机选取初始值

k=1;

m=1;

while m>delta

    
v=v0-double(subs(Fab,var,v0)/subs(dFab,var,v0)); %牛顿迭代公式

    m=norm(v-v0);

    v0=v;

    k=k+1;

    if(k>
100)

        break;%迭代步数太多，可能不收敛！

    end

end

K(i)=k;

if k<100

    t=t+1;

    V(t,:)=v; %记录收敛到的a b 值

end

end

%删除V中0元素

A=V(:,1);

B=V(:,2);

Al=logical(A);

Bl=logical(B);

VV=[A(Al) B(Bl)];

for n=1:length(VV)

   figure

   f=VV(n,1)*exp(VV(n,2)*x);

   plot(x,f,'b')

   hold on

   plot(x,y,'r')

  title(['a=',
num2str(VV(n,1)),'   ','b=',
num2str(VV(n,2))]);

  legend('y=a*e(bx)')

  xlabel('x')

  ylabel('y')

end

%均方误差

Z(1,length(VV))=0;

for j=1:length(VV)

for i=1:length(x)

   
Z(j)=Z(j)+(VV(j,1)*exp(VV(j,2)*x(i))-y(i))^2;

end

  Z(j)=sqrt(Z(j)）;

end

Zm=min(Z);%最小均方误差

 col=find(Z==Zm);

a1=VV(col,1);%使均方误差为最小的a

b1=VV(col,2);%使均方误差为最小的b

%使均方误差为最小的a，b对应的最大偏差Pd

P(length(x))=0;

for i=1:length(x)

    P(i)=abs(a1*exp(b1*x(i))-y(i));

end

Pd=max(P);

fprintf('a=%fl\n',a1)

fprintf('b=%fl\n',b1)

fprintf('均方误差:%fl\n',Zm)

fprintf('最大偏差:%fl\n',Pd)