非线性系统的Newton-Raphson和Gauss-Seidel

Question

我需要帮助弄清楚如何在Matlab中结合Newton-Raphson和Gauss-Seidel方法来求解非线性方程组。这是我的逻辑：

1. 通过求方程的导数将系统线性化并放入矩阵中。转移系统，使我们有[偏导数矩阵][H]=[原始系统]，其中H是我要求解的。H是步长
2. 首先将我对x的初始猜测插入第一个和最后一个矩阵，然后每次求解H时，我都会将该步长添加到之前对x的猜测中，以创建一个新的x。然后，我将插入更新后的x，找到另一个H，并继续该过程
3. 每次迭代都有一个新旧猜测之间的错误。当误差小于公差时，会报告这些值

我以前对Newton-Raphson和Gauss-Seidel方法进行过编码，但分别进行了编码。我的主要困惑来自于如何准确地将其设置在一起，并使用两个不同的变量。首先，我想知道我的逻辑是否正确，这个设置的pseduo代码将非常有用。

我试着对它进行了编码，我相信矩阵是正确的，我认为高斯-塞德尔函数是正确的。但我对如何正确更新猜测感到非常困惑。现在，一旦我消除了错误，它就会输出0的迭代。我被告知对线性方程做[0,0]的猜测，对非线性方程做[2,2]的猜测。基本上，我不知道如何将Newton-Raphson融入其中。我以前对NR方法进行过编码，但只使用了1个变量，所以我不明白如何将这个逻辑转换到这个代码中。这是我单独的Newton-Raphson代码：

Answer 1

您正走在正确的轨道上，但我们可以在您的代码中改进以下几点：

1. 首先，注意Newton-Raphson方法和Gauss-Seidel方法都是求解方程组的数值方法。不惜一切代价避免使用符号变量。它不仅速度慢得多，而且不能解决所有的数值问题
2. 使用和滥用MATLAB语法来处理矩阵和向量

您的问题是耦合Newton-Raphson和Gauss-Seidel解算器。我不会解决整个问题，正如你所说的，你的高斯塞德尔似乎在起作用。

我也会稍微改变你提供的一些符号，这样你就可以完全看到MATLAB易于使用的语法。

---

非线性方程组：

所提供的非线性方程组可以用向量形式写成：

function fx = F(x)
fx = zeros(2,1); % pre-allocate to create a column-vector
fx(1) = x(1)-x(2)+1;
fx(2) = x(1)^2+x(2)^2-4;
end

这创建了作为变量x1(x(和x2(y(的函数的列向量函数F。

Newton-Raphson算法需要对Jacobian进行评估。你需要计算这些导数。功能很简单，所以您可以手动执行。创建一个函数来计算作为变量函数的雅可比。

function J = jacob(x)
J = zeros(2,2);
J(1,1) = 1; % df/dx
J(1,2) = -1; % df/dy
J(2,1) = 2*x(1); % dg/dx
J(2,2) = 2*x(2); % dgdy
end

你关于牛顿-拉斐逊的逻辑是正确的。让我们用它来编码我们的Newton-Raphson:

function x = newton(fun,jacobian,x0)
% This function uses the Newton-Raphson method to find a root for the
% nonlinear system of equations F, with jacobian J, given an initial
% estimate x0.
% Here we take advantage of MATLAB syntax to easily implement our
% routine. fun and jacobian are function handles to the equations and 
% its derivatives. So it is computationally much more efficient to use
% it instead of working with symbolic variables.
% Following the steps you provided and the Newton-Raphson for one variable:
TOL = 0.00001;
MAXITER = 100;
err = 1;
iter = 0;
x = x0;
while err>TOL && iter<MAXITER % repeat until error is smaller than tolerance
% step 1: linearize the system
fx = fun(x);
J = jacobian(x);
% notice that the system is already in matrix form.
% ********************************************* %
% step 1.2: solve for stepsize H              * %
H = -Jfx; % SOLVE LINEAR SYSTEM OF EQUATIONS * %
% ********************************************* %
% step 2: add stepsize to the previous guess of x
x = x+H;
% step 3: calculate error
err = norm(H);
end
end

上述函数使用Newton-Raphson方法来求解非线性方程组的零点。让我们检查一下它是否在我们的主程序上工作：

clc; clear; close all force;
x = newton(@F,@jacob,[2;2]) % implemented Newton-Raphson
x0 = fsolve(@F,[2;2]) % Matlab built-in function
norm(x-x0) % compare both solutions

ans =
8.9902e-10

我们的函数与MATLAB的内置函数之间的差异非常小。看来这个功能是有效的。

请注意上面突出显示的步骤1.2。这是重要的一步。在这一步中，我使用MATLAB齿隙算子来求解线性系统Ax=b。以下语句具有相同的功能(求解线性方程组(：

x = AB
x = mldivide(A,B)

如果你必须使用高斯-塞德尔方法来求解线性方程组，我将把修改留给你。在步骤1.2中，你将更改高斯-塞德尔函数的齿隙算子：

H = -Jfx;
H = gaussSeidel(-J,fx);

正确地编码gaussSeidel函数(将代码转换为函数(：

function x = gaussSeidel(A,b)
% Solves the system of linear equations Ax=b by the Gauss-Seidel
% method.
%
% ...
%
end

这就行了。