matlab 非线性方程求解算法实现

非线性方程的求解算法实现

不动点法

将方程 改写成等价的形式

选择一个初始近似值 ,将它代入上式右端，即可求得

如此反复迭代计算

代码：

这次选用的函数为

改成等价形式为(我们构造迭代公式的时候，应该使迭代公式的一阶导的绝对值小于1)

%此为不动点法
syms x;
g(x)=(-7+x)^(1/5);%由于原函数为x^5-x-7=0，根据原函数求出迭代公式
exp=10^(-5);%给定判定收敛的条件
x=[];%初始化
x(1)=1.5;%给定初始点
 fprintf("迭代结果为\n%f\n",x(1));
i=1;
while(1)
    x(i+1)=g(x(i));%根据公式计算
    if(abs(x(end)-x(end-1))<exp)%如果满足收敛条件，则退出
        break;
    end
    fprintf("%f\n",x(i+1));
    i=i+1;
end

结果：

1.175596473384379 + 0.805696217928811i，一共迭代了五次

牛顿法

设已知方程有近似根,假定,将函数在点进行泰勒展开，有

于是方程可近似表示为

这是个线性方程，记其根为,则的计算公式为

这就是牛顿法公式

它的优点是收敛速度快，但是缺点是结果受初始值影响大，且每一次迭代都要算一阶导数计算麻烦，为了简便计算就有了简化牛顿法。

代码：

这次选用的函数为

clear,clc;
syms x;
%此为牛顿法
f(x)=x^2-10;%原函数
f1(x)=diff(f(x),1);%求一阶导
g(x)=x-f(x)/f1(x);%牛顿法
i=1;
esp=10^(-5);%给出判定收敛的条件
x=[];%初始化
x(1)=1;%初始点
while(1)
    x(i+1)=g(x(i));%根据迭代公式计算
    if(abs(x(end)-x(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
end

结果：

1	2	3	4	5	6	7
1	5.5	3.659090909090909	3.196005081874647	3.162455622803890	3.162277665175675	3.162277660168380

迭代了6次，最终结果为 3.16227766016838

简化牛顿法

为了简化计算，用代替,即

代码：

选用的函数同样为，来对比结果

syms x1;
%此为简化牛顿法
i=1;
f(x1)=x1^2-10;%原函数，与牛顿法的相同
f1(x1)=diff(f(x1),1);%求一阶导
h(x1)=x1-f(x1)/f1(x1(1));%简化牛顿法
x1=[];%初始化
x1(1)=1;%初始点
while(1)
    x1(i+1)=h(x1(i));%根据迭代公式计算
    if(abs(x1(end)-x1(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
end

结果：

1	2	3	4	5	6	7
1	5.5	3.65909090909091	3.19600508187465	3.16245562280389	3.16227766517568	3.16227766016838

同样迭代了六次，结果为3.16227766016838，与牛顿法的结果一样

牛顿下山法

下山法保证逐步收敛，即：

牛顿法与下山法结合起来使用，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度.

与前一步的近似值$x_k$适当加权平均作为新的改进值

其中,为下山因子

下山因子：选择下山因子时从开始，逐次减半进行试算，直到下降条件成立为止

代码：

这次选用的函数为,下面对比用下山法和不用下山法的区别，为了体现出区别，顺利进入下山法的判断条件，我的初始点选择了-9.6，离根（1.236505703391499）相去甚远。

不用下山法，只是牛顿法：

syms x2;
%此为牛顿法
f(x2)=x2^5-x2^3-1;%原函数
f1(x2)=diff(f(x2),1);%求一阶导数
y(x2)=x2-f(x2)/f1(x2);
i=1;
x2=[];%初始化
x2(1)=-9.6;%初始点
a=1;
while(1)
    x2(i+1)=y(x2(i));
 fprintf("第%d次迭代，迭代结果 %d\n",i,x2(i+1));
    if(abs(x2(end)-x2(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
end
 fprintf("\n");

结果：

迭代了25次，收敛速度比较慢，而且我们看第12次迭代，它离更偏远了，前十一次都在不断靠近，且计算和的函数值发现它并不是单调递减，所以在这里我们要用下山法改进第12次迭代的结果。

牛顿下山法结合：

代码：

syms x2;
%此为牛顿下山法
f(x2)=x2^5-x2^3-1;%原函数
f1(x2)=diff(f(x2),1);%求一阶导数
y(x2)=x2-f(x2)/f1(x2);%牛顿下山
i=1;
x2=[];%初始化
x2(1)=-9.6;%初始点
a=1;
while(1)
    x2(i+1)=y(x2(i));
    if(abs(f(x2(i+1)))>=abs(f(x2(i))))%如果后一项的函数值绝对值大于它上一个的函数值
        lamda=1;%下山因子
        while(a)%a为判断条件
            if(abs(f(y(lamda*x2(i+1)+(1-lamda)*x2(i))))<abs(f(x2(i))))%如果含有下山因子的新值对应的函数值比上一个的函数值小，说明调整成功，保持了函数值单调递减
                x2(i+1)=lamda*x2(i+1)+(1-lamda)*x2(i);%那么用含有下山因子的新值代替原本的x值
                a=0;%a=0表明退出循环
            else%不满足if条件，继续调整
                lamda=lamda/2;%lamda减半
            end
        end
    end
 fprintf("第%d次迭代，迭代结果 %d\n",i,x2(i+1));
    if(abs(x2(end)-x2(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
end
 fprintf("\n");

结果：

可以看到，第12次结果被改进，迭代结果逐步靠近，只用迭代16次就可以出结果，lambda下山因子为0.25

重根情况

知道根的重数：

不知道根的重数：

代码：

本次选用的函数为,并不知道根的重数，所以套用第二条公式。初始值为1.5

syms x3;
%重根情况，不知道几重
f(x3)=x3^6-3*x3^4+1;%原函数
f1(x3)=diff(f(x3),1);%求一阶导数
f2(x3)=diff(f(x3),2);%求二阶导数
k(x3)=x3-(f(x3)*f1(x3)/((f1(x3))^2-f(x3)*f2(x3)));%带入公式
 x3=[];%初始化
 i=1;
 x3(1)=1.5;%初始点
 while(1)
    x3(i+1)=k(x3(i));
    fprintf("第%d次迭代，迭代结果 %d\n",i,x3(i+1));
    if(abs(x3(end)-x3(end-1))<esp)
        break;
    end
    i=i+1;
 end
 fprintf("\n");

结果：

弦截法

特点：两个初始值，一阶导数用差商替代

收敛阶：

代码：

本次选用的函数为，初始值为0.7和0.6。

syms x4;
%此为弦截法
f(x4)=(x4^2)*exp(x4)-10;
 x4=[];
 i=2;
 x4(1)=0.7;%给定两个初始点
 x4(2)=0.6;
 while(1)
    a=x4(i)-f(x4(i))*(x4(i)-x4(i-1))/(f(x4(i))-f(x4(i-1)));%代入公式
    x4(i+1)=a;
     fprintf("第%d次迭代，迭代结果 %d\n",i-1,x4(i+1));
    if(abs(x4(end)-x4(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
 end
 fprintf("\n");

结果：

迭代了14次

抛物线法

特点：三个初始值，运用抛物线的交点作为根的近似位置

收敛阶：

代码：

syms x5;
 %此为抛物线法
f(x5)=(x5^2)*exp(x5)-10;
 x5=[];
 i=3;  
 x5(1)=0.7;
 x5(2)=0.6;%给定初始点
 x5(3)=0.56532;
 while(1)
     f_yijie_21=(f(x5(i-1))-f(x5(i-2)))/(x5(i-1)-x5(i-2));%求一阶差商f[x2,x2]
     %fprintf("一阶%f\n",f_yijie_21);
      f_yijie_31=(f(x5(i-2))-f(x5(i)))/(x5(i-2)-x5(i));%求一阶差商f[x3,x1]
      %fprintf("一阶%f\n",f_yijie_31);
      f_yijie_32=(f(x5(i-1))-f(x5(i)))/(x5(i-1)-x5(i));%求一阶差商f[x3,x2]
     % fprintf("一阶%f\n",f_yijie_32);
     f_erjie=(f_yijie_31-f_yijie_32)/(x5(i-2)-x5(i-1));%求二阶差商f[x2,x1,x0]
     %fprintf("二阶%f\n",f_erjie);
     omega=f_yijie_32+f_erjie*(x5(i)-x5(i-1));% 求w
     %fprintf("omega%f\n",omega);
     if omega>0%判断omega符号
          a=x5(i)-((2*f(x5(i)))/(omega+(omega^2-4*f(x5(i))*f_erjie)^(1/2)));%计算结果
     else
         a=x5(i)-((2*f(x5(i)))/(omega-(omega^2-4*f(x5(i))*f_erjie)^(1/2)));
     end
    x5(i+1)=a;
     fprintf("第%d次迭代，迭代结果 %d\n",i-2,x5(i+1));
    if(abs(x5(end)-x5(end-1))<esp)%如果满足收敛条件，则退出
        break;
    end
    i=i+1;
 end

结果：

迭代了6次，迭代次数比弦截法的少，迭代效果比弦截法好，与理论相符。