使用条件
二分法需要知道两个自变量,分别是一个根的两侧
牛顿法迭代法是最常用的方法,收敛性信赖于初值,取不同的初值可以的方程不同的根,函数用的是一阶导数,输入的是一个猜想的可能的值
割线法给定两个初值再带入计算,比如要在2附近求一个根,那就可以假设这个范围是(1.9,2)
拟牛顿法这个比较方便,用时最好可以找到一个好的初始值
Halley法需要知道函数值以及它的一阶求导、二阶求导
这里我就从计算代码的角度来讲解,在下面也会按照上面这个顺序给出代码,遇到方程组直接带入已知条件就可以得到答案。
二分法
基本函数是这样子的:y = dichotomy(fun,a,b,tol);二分法的算法要输入四个变量,fun,a,b,tol:函数,一个根的左右点,tol=1.0e-6
function y = fun(x)
y = x^3 - 5 * x + 4.272;
上面这个就是定义的fun,每次的输入的方程不同,第一条不动,直接改第二行就可以的。比如这里我们要计算的方程y = x^3 - 5 * x + 4.272;
我们是可以通过简单计算得到一个根的两侧分别是1和1.3
那在窗口指令指令中输入x=dichotomy(’fun‘,1,1.3,1.0e-6)就可以得到结果
function y = dichotomy(fun,a,b,tol)
if nargin < 4
tol = 1.0e-5;
end
n = 1;
if feval(fun,a)*feval(fun,b)<0
c = (a+b)/2;
while (abs(b-c)>tol) && (abs(feval(fun,c))>tol)
if (feval(fun,c)*feval(fun,a)>0)
a = c;
c = (a+b)/2;
elseif (feval(fun,c)*feval(fun,a)<0)
b = c;
c = (a+b)/2;
else
y = c;
tol = 100;
end
n = n + 1;
end
y = c;
elseif feval(fun,a)==0
y = a;
elseif feval(fun,b)==0
y = b;
else
disp('there may not be a root in the interval');
end
n
function y = fun(x)
y = x^3 - 5 * x + 4.272;
牛顿法
还是用刚才那道题,y = x^3 - 5 x + 4.272,一阶导是y = 3 x^2 - 5;
function y = dfun(x)
y = 3 * x^2 - 5;
下面的是具体的算法,根据x = newton(x0,tol),我们只需要输入一个我们猜想的值就可以。但是有一定的误差
function x = newton(x0,tol)
if nargin < 2
tol = 1.0e-5;
end
x = x0 - fun(x0)/dfun(x0);
n = 1;
while (norm(x-x0)>tol) && (n<1000)
x0 = x;
x = x0 - fun(x0)/dfun(x0);
n = n + 1;
end
n
割线法
这里我们用割线法求y = x^3 - 5 * x + 4.272在方程x=2的根,输入上要用两个初始值,比如说现在来计算就可以输入x=secant(2,1.9,10e-6)
function x = secant(x0,x1,tol)
if nargin < 3
tol = 1.0e-5;
end
x = x0 - fun(x0) * (x0 - x1) / (fun(x0) - fun(x1));
n = 1;
while (abs(x0-x1) > tol) && (n <= 1000)
x1 = x0;
x0 = x;
x = x0 - fun(x0) * (x0 - x1) / (fun(x0) - fun(x1));
n = n + 1;
end
n
拟牛顿法
这里我们可以直接找到一个初始值输入,比如说broyden2,10e-6),如果不知道不确定也没关系,至少要知道一个范围。比如说给个范围(0.5,0.5)有下面这个函数
function y = funm(x)
y(1,1) = x(1,1) - 0.7 * sin(x(1,1)) - 0.2 * cos(x(2,1));
y(2,1) = x(2,1) - 0.7 * cos(x(1,1)) + 0.2 * sin(x(2,1));
那就可以输入x = broyden(x0,tol)
function x = broyden([0.5,0.5],tol)
if nargin < 2
tol = 1.0e-5;
end
A = eye(size(x0,1));
x = x0 - A \ funm(x0);
n = 1;
while (norm(x - x0) > tol) && (n < 1000)
x0 = x;
x = x0 - A \ funm(x0);
p = x - x0;
q = funm(x) - funm(x0);
A = A + (q - A*p)*p'/norm(p)^2;
n = n + 1;
end
n
Halley法
这个要求二阶导,比如说第一个道题,y = x^3 - 5 * x + 4.272;,二阶导数是下面这个输入
function y = d2fun(x)
y = 6 * x;
最后输入x = halley(1,10e-6)就可以计算出一个结果
function x = halley(x0,tol)
if nargin < 2
tol = 1.0e-5;
end
m = size(x0,1);
x = x0 - (eye(m) - 1/2 * (dfun(x0) \ d2fun(x0)) * (dfun(x0) \ fun(x0))) \ ...
(dfun(x0) \ fun(x0));
n = 1;
while (norm(x - x0) > tol) && (n < 1000)
x0 = x;
x = x0 - (eye(m) - 1/2 * (dfun(x0) \ d2fun(x0)) * (dfun(x0) \ fun(x0))) \ ...
(dfun(x0) \ fun(x0));
n = n + 1;
end
n