MATLAB求解非线性方程组（牛顿拉夫逊方法）

Equation.m 函数

这个函数用来写非线性方程组，注意整理成标准形式

function [f,variable]=Equation(~) 
%定义非线性方程组,方程组右端项为0  
syms x1 x2  x3 x4 %x5 x6 x7 x8 x9 x10 x11 x12
f1=x1+x1*x2+x3-8*x4-8;
f2=x3+x2+x1*x3+x4*x1+x4;
f3=x1*x1+x3*x4-5+x3+x4;
f4=x2+x4+x3*x4;
%f5=x6*x5*x4;
%f6=

%f=[f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 ];
f=[f1 f2 f3 f4];
variable=[x1 x2 x3 x4];
end

牛顿拉夫逊方法迭代求解

clear
[f,variable]=Equation()

N=5;%迭代次数
[row,size]=size(variable);
x0=zeros(size,1); %初值
for i =1:N
%     delta_y=zeros(size,1)-double(subs(Equation(),...
%     {'x1' 'x2' 'x3' 'x4' 'x5' 'x6' 'x7' 'x8' 'x9' 'x10' 'x11' 'x12'},...
%     {x0(1) x0(2) x0(3) x0(4) x0(5) x0(6) x0(7) x0(8) x0(9) x0(10) x0(11) x0(12)})); 
%     J=double(subs(jacobian(f,variable);...
%       {'x1' 'x2' 'x3' 'x4' 'x5' 'x6' 'x7' 'x8' 'x9' 'x10' 'x11' 'x12'},...
%     {x0(1) x0(2) x0(3) x0(4) x0(5) x0(6) x0(7) x0(8) x0(9) x0(10) x0(11) x0(12)}));  
    delta_y=zeros(size,1)-double(subs(Equation(), ...
    {'x1' 'x2' 'x3' 'x4' }, ...
    {x0(1) x0(2) x0(3) x0(4)}))'; 
    J=double(subs(jacobian(f,variable), ...
      {'x1' 'x2' 'x3' 'x4' }, ...
    {x0(1) x0(2) x0(3) x0(4)}));
    delta_x=J\delta_y
    x0=x0+delta_x;
    
    if(norm(delta_y,"inf")<0.5)%判定收敛
        disp(i);
        disp(x0);
        break;
    end
end

本例运行结果输出：

[x1 + x3 - 8*x4 + x1*x2 - 8, x2 + x3 + x4 + x1*x3 + x1*x4, x1^2 + x3 + x4 + x3*x4 - 5, x2 + x4 + x3*x4]
[x1, x2, x3, x4]
delta_x = 4×1    
    48
    -5
     0
     5

delta_x = 4×1    
  -23.9728
    2.5888
   -0.0210
   -2.4835

delta_x = 4×1    
  -11.9224
    1.5449
   -0.2289
   -1.0432

delta_x = 4×1    
   -4.9237
   10.1295
  -12.8115
   11.3410

delta_x = 4×1    
   -1.9197
    1.8696
    5.1976
   -6.3692

delta_x = 4×1    
   -0.7955
   -2.1395
    2.5462
   -2.7440

delta_x = 4×1    
   -0.3771
   -1.5859
    1.1497
   -0.9999

delta_x = 4×1    
   -0.1128
   -0.4665
    0.3133
   -0.2430

delta_x = 4×1    
   -0.0095
   -0.0376
    0.0249
   -0.0184

     9
    3.9665
    6.9032
   -3.8297
    2.4397

Matlab 算是一个比较好的计算工具了，偶尔可以用用，现解方程组如下：

syms a b c
f1 = a*exp(0)-c
f2 = a*exp(481*b)-c-2107/16
f3 = a*exp(962*b)-c-2107
[a b c] = solve([f1 f2 f3])

之后输出：

a =
2107/(exp(2log(15)) - 1)
b =
log(15)/481
c =
2107/(exp(2log(15)) - 1)

ok 啦！

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;

上面这个就是定义的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