MATLAB求解方程的根

%% 2022-1-12求方程的根
%{
求解方程的根 x^3+2^x+1=0，分别采用二分法，牛顿法，简单迭代法，割线法，
上坡法和下坡法以及zeroin算法
%}
%% x^3+2^x+1=0 判断方程的根数，以及根的范围
clc;close all;clear all;
%画函数图像判断
syms ('f','x');
f=@(x)x.^3+2.^x+1;%就是一个表达式就这样写了，方便
h=figure('name','函数');
fplot(f,'k');
magnifyOnFigure(h, 'displayLinkStyle','straight',...
                    'EdgeColor', 'black',...
                    'secondaryAxesYLim',[-5,5],...
                    'secondaryAxesXLim',[-3,3],...
                    'magnifierShape', 'rectangle',...
                     'mode', 'interactive',...
                     'initialPositionMagnifier',[200 180 180 40],...
                    'frozenZoomAspectratio', 'off',...
                     'edgeColor','g',...
                    'edgeWidth', 2);
%通过图像的局部放大图可以，发现根的范围在[-3,3]

%判断函数有几个根
figure('name','一阶导数');
diff(f,x,1);%log函数是以e为底
df=@(x)(2.^x.*log(2) + 3.*x.^2);%这里不加点会报错
fplot(df,[-3,3]);
%通过导函数的图像可以发现这个函数是单调递增的，根据零点定量+单调可知，只有一个根
%% 1.通过二分法求解，土木一般精度没有那种特别高，我们让它有四位有效数
epsilon=0.5*10^(-3);b=3;a=-3;k=0;
while (b-a)/(2^(k+1))>epsilon
    k=k+1;
end
%f(a)<0,f(b)>0
lef=a;rig=b;
for i=1:k
    middle=(rig+lef)/2;
    if f(middle)*f(rig)>0
        rig=middle;
    elseif f(middle)*f(lef)>0
         lef=middle;
    end
end
disp(strcat('二分法方程的根为：',num2str(middle)));
%% 2.通过简单迭代法求解
% x=-(2^x+1)^(-3)
h=@(x)(-(2.^x+1).^(1/3));
diff(h,x,1);
Dh=@(x)-(2.^x.*log(2))./(3.*(2.^x + 1).^(2/3));
figure(3);
fplot(Dh)
%通过图像可以看出h(x)的导数在0-(-1)之间,h(a),h(b)都在[a,b]内，故迭代格式收敛
xj=2; %任取一个初始点
xi=h(xj);
while abs(xi-xj)>epsilon
    xj=xi;
    xi=h(xj);
end
disp(strcat('简单迭代法方程的根为：',num2str(xi)));
%{
可以发现二分法求得的解为-1.1331，简单迭代法求得的解为-1.1334，
这和四位有效数的要求是满足的
%}

%% 牛顿迭代法大范围不收敛，所以用局部收敛性
figure('name','二阶导数');
diff(f,x,2);
df=@(x)(6.*x + 2.^x*log(2).^2);%这里不加点会报错
fplot(df,[-3,3]);
xp=-1.133400;
xq=xp-(f(xp)/df(xp));
while abs(xq-xp)>epsilon*10^(-1)
 xp=xq;
 xq=xp-(f(xp)/df(xp));
end
disp(strcat('牛顿迭代法局部方程的根为：',num2str(xq)));

%% 割线法
x1=-1.13331;
x2=-1.13334;
dif=(f(x2)-f(x1))/(x2-x1);
x3=x2-(f(x2)/dif);
while  abs(x3-x1)>epsilon*10^(-1)
 x1=x2;
 x2=x3;
 dif=(f(x2)-f(x1))/(x2-x1);
 x3=x2-(f(x2)/dif);
end
disp(strcat('割线法局部方程的根为：',num2str(x3)));
%% steffenson这边称为上坡法
%% 上坡法
xd=-1.13331;
difd=(f(xd+f(xd))-f(xd))/f(xd);
xu=xd-(f(xd)/difd);
while  abs(xu-xd)>epsilon*10^(-1)
 xd=xu;
 difd=(f(xd+f(xd))-f(xd))/f(xd);
 xu=xd-(f(xd)/difd);
end
disp(strcat('上坡法局部方程的根为：',num2str(xu)));
%% MATLAB一些自带的求根函数

%% Zerions算法将二分法的可靠性和割线法快速性结合
options = optimset('Display','final','FunValCheck','off','TolX',epsilon);
z=fzero(f,[-3,3],options);
disp(strcat('Zerions算法方程的根为：',num2str(z)))
%% 寻找单变量的局部极小值，在求解根的时候也是非常有用的这里，演示一下
[x,fval]=fminbnd(f,-3,3,options);
%fminsearch使用无导数法计算无约束的多变量函数的最小值
%这里主要应用的是使用Lagarias 等的单纯形搜索法。
%% 最后看一下各个算法的结果和误差
A={};
A{1,1}='二分法';A{1,2}=middle;A{1,3}=f(middle);
A{2,1}='简单迭代法';A{2,2}=xi;A{2,3}=f(xi);
A{3,1}='牛顿迭代法局部';A{3,2}=xq;A{3,3}=f(xq);
A{4,1}='割线法局部';A{4,2}=x3;A{4,3}=f(x3);
A{5,1}='牛顿上坡法局部';A{5,2}=xu;A{5,3}=f(xu);
A{6,1}='Zerions算法';A{6,2}=z;A{6,3}=f(z);
Zx=A(:,1);
n=size(A);B=[];
for i=1:n(1)
    B(i,1)=abs(A{i,3}-0);
end
figure
plot([1:2:12],B,'-*k','MarkerEdgeColor','r','MarkerSize',5);
set(gca,'XTick',1:2:12,'XTickLabel',Zx,'XTickLabelRotation',45);
title('各方法误差分析图')