非线性方程组数值求解算法——MATLAB源码

---

目录

关于什么？

1、牛顿法

2、拟牛顿法

---

---

关于什么？

主要是牛顿法和拟牛顿法求解非线性方程组的算法步骤以及MATLAB代码实现。编程的思路是主要的，CV大法虽然快，但是重要的还是仔细研究算法步骤以及实现过程。

以二变量的二次非线性方程组求解为例，进行示例。有问题在评论区留下，我来解答。

推荐参考书籍：

[1] 薛毅. 数值分析与科学计算[M]. 科学出版社, 2011.

---

1、牛顿法

开门见山，原理性的东西不做解释，可以查阅书本，直接3，2，1，上算法步骤！

参考书例5.2：

---

MATLAB 编程实现（初次接触建议细看，当然编程方式不唯一）：

%% 非线性方程组 - Newton 解法
% Create Time : 2022 04 17
% Finish Time : 2022 04 17
% Author      : Xu Y. B. (CSDN ID:在路上，正出发)

% Ref         :[1] 薛毅. 数值分析与科学计算[M]. 科学出版社, 2011. 5.2节

%% CLEAR
clc;
clearvars;
close all;

%% 符号及方程组定义
syms x1 x2;
f1(x1,x2) = x1^2+x2^2-5;
f2(x1,x2) = (x1+1)*x2-(3*x1+1);

F(x1,x2) = [f1;f2];

f1_x1_Pd(x1,x2) = diff(f1,x1);% f1 对 x1 求偏导
f1_x2_Pd(x1,x2) = diff(f1,x2);% f1 对 x2 求偏导
f2_x1_Pd(x1,x2) = diff(f2,x1);% f2 对 x1 求偏导
f2_x2_Pd(x1,x2) = diff(f2,x2);% f2 对 x2 求偏导

J(x1,x2)=[f1_x1_Pd,f1_x2_Pd;...
          f2_x1_Pd,f2_x2_Pd];% F(x)的 Jacobi 矩阵
     
%% 迭代求解
% 参数设置
N = 100;% 设置最大迭代次数
epsilon = 1e-3;% 设置迭代的精度
x_start = [1;1];% 迭代初始点
x_k_solve = zeros(2,N+1);% 存放每次迭代结果，以便分析，第一列为初值
x_k_solve(:,1) = x_start;
dk = zeros(2,N);% 存放每次迭代的 dk
for i = 1:N
    dk(:,i) = J(x_k_solve(1,i),x_k_solve(2,i))\(-F(x_k_solve(1,i),x_k_solve(2,i)));
    if(sqrt(sum(dk(:,i).^2)) < epsilon)
        x_k_solve(:,i+1) = x_k_solve(:,i) + dk(:,i);
        break;
    else
        x_k_solve(:,i+1) = x_k_solve(:,i) + dk(:,i);
    end
end     

计算结果：

 课本结果：

 二者完全一致，验证了程序的正确性。

2、拟牛顿法

算法步骤：

参考书的例题5.3

---

 MATLAB代码：

%% 非线性方程组 - 拟 Newton 解法
% Create Time : 2022 04 17
% Finish Time : 2022 04 17
% Author      : Xu Y. B. (CSDN ID:在路上，正出发)


% Ref         :[1] 薛毅. 数值分析与科学计算[M]. 科学出版社, 2011. 5.2节

%% CLEAR
clc;
clearvars;
close all;

%% 符号及方程组定义
syms x1 x2;
f1(x1,x2) = x1^2+x2^2-5;
f2(x1,x2) = (x1+1)*x2-(3*x1+1);

F(x1,x2) = [f1;f2];

     
%% 迭代求解
% 参数设置
N = 100;% 设置最大迭代次数
A0 = [1,0;0,1];% 初始矩阵
epsilon = 1e-3;% 设置迭代的精度
x_start = [1;1];% 迭代初始点
x_k_solve = zeros(2,N+1);% 存放每次迭代结果，以便分析，第一列为初值
x_k_solve(:,1) = x_start;
sk = zeros(2,N);% 存放每次迭代的 sk
Ak = A0;% 迭代矩阵初始化
y = zeros(2,N);% 矩阵y初始化
for i = 1:N
    sk(:,i) = Ak\(-F(x_k_solve(1,i),x_k_solve(2,i)));
    if(sqrt(sum(sk(:,i).^2)) < epsilon)
        x_k_solve(:,i+1) = x_k_solve(:,i) + sk(:,i);
        break;
    else
        x_k_solve(:,i+1) = x_k_solve(:,i) + sk(:,i);
        y(:,i) = F(x_k_solve(1,i+1),x_k_solve(2,i+1))-F(x_k_solve(1,i),x_k_solve(2,i));
        Ak = Ak+((y(:,i)-Ak*sk(:,i))*(sk(:,i).'))/((sk(:,i).')*sk(:,i));
    end
end 

运算结果：

 课本结果：

 由结果可以明显看出，二者计算结果一致，证明程序的准确性。