MATLAB求解非线性方程组的五种方法

MATLAB求解非线性方程组的五种方法

求解线性方程分为两种方法–二分法和迭代法
常见的方法一共有5种
二分法
迭代法
牛顿法
割线法
拟牛顿法
Halley法

使用条件
二分法需要知道两个自变量，分别是一个根的两侧

牛顿法迭代法是最常用的方法，收敛性信赖于初值，取不同的初值可以的方程不同的根，函数用的是一阶导数，输入的是一个猜想的可能的值

割线法给定两个初值再带入计算,比如要在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;
1
2

上面这个就是定义的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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

function y = fun(x)
y = x^3 - 5 * x + 4.272;
1
2

牛顿法
还是用刚才那道题，y = x^3 - 5 * x + 4.272，一阶导是y = 3 * x^2 - 5;

function y = dfun(x)
y = 3 * x^2 - 5;
1
2

下面的是具体的算法，根据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
1
2
3
4
5
6
7
8
9
10
11
12

割线法

这里我们用割线法求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
1
2
3
4
5
6
7
8
9
10
11
12
13

拟牛顿法

这里我们可以直接找到一个初始值输入，比如说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));
1
2
3

那就可以输入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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Halley法

这个要求二阶导，比如说第一个道题，y = x^3 - 5 * x + 4.272;，二阶导数是下面这个输入

function y = d2fun(x)
y = 6 * x;
1
2

最后输入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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16