MATLAB 数值求解方程根

本文介绍如何利用MATLAB求解给定区间内方程的数值解

1. 问题

针对如下函数，求[0,1]区间内的数值解（该方程无解析解）
$$f(x)=-(-8x^7+3x^2)\sin (3\pi (x-1))-3\pi (-8x^8+x^3)\cos (3\pi (x-1))=0$$
通过绘图可以看出，该方程有5个根，除端点外，其余三个根分别落于[0,0.4],[0.4,0.6]和[0.6,1]内，注意MATLAB内置的fsolve函数只能处理单调区间内的零点问题，因此需要先自动搜寻到这三个区间，所依据的原理即零点定理。

图 1 示例函数

2. 解答

思路：先利用符号函数和差分运算确定变号区间，然后分区间利用内置函数fsolve求解即可，最后需要对端点进行检查。
本案例的解为：[0,0.2605,0.5514,0.8300,1.0000]
MATLAB 代码如下：

clc;
clear all;
close all;
solution_range = [0,1]; % 定义域
%% 1. 可视化
a = linspace(solution_range(1),solution_range(2),100);
y = myfun(a);
plot(a,y,'k-')
hold on
plot(a,0*y,'r-.')
xlabel('x')
ylabel('y')
title('函数图像')
%% 2.求解
% 2.1 零点搜索
[~,pos] = find(abs(diff(heaviside(y)))); % 寻找零点所在位置
xdomain = []; % 初始化
for i = 1:length(pos)
    xdomain = [xdomain;a(pos(i))]; % 更新
end
% 2.2 正式求解
x = [];
for j = 1:length(xdomain)-1
    x_latent = fsolve(@myfun,[xdomain(j),xdomain(j+1)]); % 数值解
    x = [x;x_latent(end)]; % 更新数值解集合
end
% 2.3 端点检查
if ~ismember(x,solution_range(1))
    x = [solution_range(1);x]; % 加入左端点
end
if ~ismember(x,solution_range(2))
    x = [x;solution_range(2)]; % 加入右端点
end
x
%% 3. 自定义函数
function f = myfun(x)
f = -(- 8.*x.^7 + 3.*x.^2).*sin(3.*pi.*(x - 1)) - 3.*(- x.^8 + x.^3).*pi.*cos(3.*pi.*(x - 1));
end

3. 参考资料

matlab函数零点求法,Matlab之函数零点 | 学步园