matlab 局部寻根,MATLAB算法の二分法全局寻根

今天是二分法介绍的最后一个版本，在《MATLAB算法の二分法改进版》中提出通过matlab矢量化计算实现类并行二分法快速寻根。有读者指出二分法是不是只能用于函数单调的区间，对于非单调区间能不能使用二分法解决问题？这里需要说明的是区间单调并不是二分法寻根的充分必要条件，只要需要在根的左右两侧的函数值是异号即可。至于在某个区间存在多根的情况，可以采用“分而治之”的办法完美解决。

先给出二分法的定义：对于在区间[a,b]上连续且单调的函数f(x)，若满足条件f(a)*f(b)<0，则函数f(x)在此区间上必存在根。通过不断把此区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点准确值或近似值的方法。

现在具体说明如何实现全局寻根，如果函数f(x)在区间[a,b]存在多个根，则在区间[a,b]上函数值f(x)必然是正负交替出现的一系列值，而正负号发生转变的位置区间内必然存在函数f(x)的根。通过这种方法就可以锁定各个根所在的大致区域，然后再使用二分法进行精细寻找。

对于寻根部分，使用《二分法改进版》中的方法来进行n分寻根，即将区间[a,b]上n等分，计算出各个等分点的函数值，将函数值正负值发生改变处的两个等分点作为新的区间[na,nb],继续将此区间n等分，重复上面步骤，将区间差绝对值Δ=|na-nb|作为精度控制标识，达到精度后，取na与nb的平均值作为函数的根。

具体实现步骤：

1、将区间[a,b]进行n等分，计算各点的函数值。若函数值全为正或全为负，则此区间不存在根；若函数值有正有负，则找出一系列正负号发生转变的位置。

2、计算根的个数，对各个根所在区域进行分别进行寻根。

2.1、n等分区间[a,b],并计算各等分点的函数值。

2.2、查找函数值中是否存在等于的点，如存在则该点即为所寻之根，若不存在则找出函数值正负发生变化的两个位置，将前者作为a，后者作为b。

2.3、通过比较区间差绝对值与计算精度e的大小来判断是否达到预设条件，若|a-b|

3、得到在区间[a,b]满足计算精度的函数f(x)的所有根。

二分法全局寻根matlab源代码

问题定义：求函数f(x) =

3*exp(sin(x.^3))-9*cos(x.^2)-5.6*sin(x)在区间[-2,3]的根，计算精度为10^-9.

clc;clear;close all;

format long;

x = -2:0.01:3;

% 定义参考y值，即y=0的直线

ty = zeros(1,length(x));

% 定义在区间[-2,3]上的连续函数fun

fun = @(x) 3*exp(sin(x.^3))-9*cos(x.^2)-5.6*sin(x);

y = fun(x);

plot(x,y,'b.-',x,ty,'r--','LineWidth',3.5);

xlabel('x轴');

ylabel('y轴');

title('二分法全局寻根测试');

locM = find(y<0);

locP = find(y>0);

if not(isempty(locM)) && not(isempty(locM))

% 定义计算精度

% 区间划分份数

divQJ = 1000;

% 查找正负号变化的位置

b=diff(y>0);

locf = find(b~=0);

locb = find(b~=0)+1;

lenL = length(locf);

% 用qJ的存储根所在区域位置

qJ(1:2:lenL*2) = x(locf);

qJ(2:2:lenL*2) = x(locb);

JG = zeros(1,lenL);

% 循环计算得出各个根

for k = 1:lenL

% 调用改进版二分法子程序

JG(k) = dichotomy(fun,qJ(2*k-1),qJ(2*k),divQJ,ep);

disp('此区间不存在根！！！');

nY = fun(JG);

% 绘制近似交点

plot(JG,nY,'r.','MarkerSize',50);

function re = dichotomy(fun,a,b,divQJ,ep)

while(tmp>ep)

tX = linspace(a,b,divQJ);

tY = fun(tX);

loc = find(tY == 0);

if isempty(loc)

locM = find(tY<0);

locP = find(tY>0);

if tY(1)<0

a = tX(locM(end));

b = tX(locP(1));

tmp = abs(a-b);

a = tX(locP(end));

b = tX(locM(1));

tmp = abs(a-b);

mid = (a+b)/2;

mid = tX(loc);

如有未尽之处，烦请大家留言指正，先行感谢！！！