MATLAB-非线性方程的数值解法——二分法

Posted 2022-11-22 Kaerou

二分法在用计算机求非线性方程解的数值方法中是最简单的一种，用人工算效率很低，但用计算机运算时还是一种很有效的方法。本文主要参考《计算方法》李大美 李素贞 朱方生编著

目录

1. 原理
2. 计算步骤
3. 程序框图
4. MATLAB实现
   4.1.按照程序框图进行编写
   4.2.先估算二分次数再进行二分
5. 例题

原理

二分法的数学理论基础是闭区间上连续函数的一个基本性质，即设f(x)在闭区间[a,b]上连续且f(a)(b)<0，则在区间内至少存在一个点α，使得f(α)=0

​​​
记a₀=a, b₀=b, 称区间[a₀, b₀]为方程f(x)=0的有根区间

对分区间[a₀, b₀]可得中点​​​​​​ x₀ 并计算出f(x₀)

若恰好有f(x₀)=0，则方程的根为


若f(x₀)≠0，则计算乘积f(x₀)f(b₀)

1.若f(x₀)f(b₀)>0，则根在区间[a₀, x₀]内，记a₁=a₀，b₁=x₀

2.若f(x₀)f(b₀)<0，则根在区间[x₀, b₀]内，记a₁=x₀，b₁=b₀

区间[a₁, b₁]是包含根的新区间，在旧的有根区间内
​

​
记d₁为区间[a₁, b₁]的长度，则

再将区间[a₁, b₁]对分，重复上述过程，可得方程新的有根区间[a₂, b₂]，其长度为

如此重复n次，若还没有找到方程的根，就得到了包含方程根的区间的一个序列：[a₀, b₀], [a₁, b₁], ···, [a_n, b_n], ···

这个闭区间序列具有这样的特点：后一个区间落在前一个区间内，且其长度只有前一个区间长度的一半

记区间[a_n, b_n]的长度为d_n，则该区间与最初区间[a, b]的长度关系为

但n充分大时，可取区间[a_n, b_n]的中点作为方程f(x)=0的一个实根α的近似值，且它们满足关系式

方程根的近似值x*的绝对误差小于最初区间长度的2ⁿ⁺¹分之一

对于给定的精度ε，可估算出二分法所需的二分次数

计算步骤

程序框图

MATLAB实现

按照程序框图进行编写

有两个允许误差ε₁，ε₂，且未估算二分次数n

clear;clc;

%使用二分法求非线性连续函数的根，无法求复数根和偶数重根
%初始未计算二分次数

a0=1;b0=2; %输入有根区间

accuracy_1=10^(-4);accuracy_2=5*10^(-4); %设定两个允许误差

y1=a0^6-a0-1; y2=b0^6-b0-1; %设定方程，若修改还需修改 20 line 处的方程

a=zeros(100,1);b=a;x=a;form=zeros(100,5); %提前创建矩阵，提高运算速度，因未计算二分次数，所以设定了100次，可适当修改
if y1*y2>0
    fprintf ('区间内有偶数个根或者无实根')
else
    n=0; %n为2分次数
    a(1)=a0;b(1)=b0;
    while ((b(n+1)-a(n+1))>accuracy_2) %循环至满足误差要求
        x(n+1)=(a(n+1)+b(n+1))/2;%求当前区间中点
        y=(x(n+1))^6-x(n+1)-1;yb=(b(n+1))^6-b(n+1)-1; %求此中点、右端点对应的f(x)
        if abs(y)<accuracy_1 %若满足第一个允许误差，即|f(x)|<ε1，填充表格并中断循环
            form(n+1,1)=n;form(n+1,2)=a(n+1);form(n+1,3)=b(n+1);form(n+1,4)=x(n+1);form(n+1,5)=y;
            break;
        else
            if y*yb>0 %y*yb大于0则让右端点等于中点，左端点保持不变，若y*yb<0则相反
                b(n+2)=x(n+1);
                a(n+2)=a(n+1);
            else
                a(n+2)=x(n+1);
                b(n+2)=b(n+1);
            end
        end
        form(n+1,1)=n;form(n+1,2)=a(n+1);form(n+1,3)=b(n+1);form(n+1,4)=x(n+1);form(n+1,5)=y; %填充表格
        n=n+1; %累计二分次数
    end
end

先估算二分次数再进行二分

先用误差估计式进行估算，求出n，再进行二分

clear;clc;

%使用二分法求非线性连续函数的根，无法求复数根和偶数重根
%初始根据误差估计式计算二分次数n

a0=1;b0=2; %输入有根区间

accuracy=5*10^(-4); %设定允许误差

y1=a0^6-a0-1;y2=b0^6-b0-1;%设定方程，若修改还需修改 20 line 处的方程

n=ceil(log2((b0-a0)/accuracy/2)); %由误差估计式计算二分次数，ceil：向正无穷取整
a=zeros(n,1);b=a;x=a;form=zeros(n,5); %根据二分次数提前创建矩阵，提高运算速度
if y1*y2>0
    fprintf ('区间内有偶数个根或者无实根')
else
    a(1)=a0;b(1)=b0;
    for i=1:n+1
        x(i)=(a(i)+b(i))/2;%求当前区间中点
        y=(x(i))^6-x(i)-1;yb=(b(i))^6-b(i)-1; %求此中点、右端点对应的f(x)
        if y*yb>0 %y*yb大于0则让右端点等于中点，左端点保持不变，若y*yb<0则相反
            b(i+1)=x(i);
            a(i+1)=a(i);
        else
            a(i+1)=x(i);
            b(i+1)=b(i);
        end
        form(i,1)=i-1;form(i,2)=a(i);form(i,3)=b(i);form(i,4)=x(i);form(i,5)=y; %填充表格
    end
end

运算结果如下图

方法一

方法二

例题

求x₃-x-1=0在[1, 1.5]内的解，要求ε<0.01

将a=1，b=1.5，accuracy=0.01，代码里的方程也改为x₃-x-1，然后运算代码，可得图表


运行所得结果与例题参考答案高度重合，因此所求方程的近似解为x*=1.3203，满足允许的误差精度