如何在 Matlab 中使用雅可比和牛顿法求解（非线性）方程组

Question

我正在尝试构建一个函数，该函数可以使用牛顿法在 Matlab 中求解具有 n 个未知数的 n-1（非线性）方程组。

上网查了一下，得出以下程序：

使用 matlab 的 'null' 函数应用于一阶偏导数 (Jacobian) 数组以获得 m 元素向量，（称此向量为 v，）它与所有 m-1 正交偏音的向量。沿着这个向量前进 h0 距离后，再次使用 'null' 函数，这次是从那个前进的点找到与 v 正交的 m-1 维子空间。将所有原始坐标投影到该子空间中的坐标，您可以执行 m-1 维 Newton-Raphson 过程，只要您的步长 h0 不太大，该过程应该会非常迅速地收敛。然后将其转换回原始坐标。

我有以下方程组：

fun = @(x,y,z)[x.^3+y.^2+z.^2,x.^2-y.^3+sin(z)]

它是雅可比行列式：

Jac = @(x,y,z)[3*x^2,    2*y,    2*z; 2*x, -3*y^2, cos(z)]

初步猜测： x0 = [1,1,1]

在我的函数中，我将矢量转换为数字，以便我可以使用它们： x = num2cell(x0)

那我要进行N次迭代：

for numIterations = 1:N

%calculate the function values at the current iteration
F = fun(x{:})

%calculate the jacobian matrix
J = Jac(x{:})

%find direction in which the curve extends:
v = null(J)

xnew = xold + h0 * v
xnew = xnew'

subspace = null(xnew)

nextX = multinewton(f,df,xnew)

现在卡在牛顿法计算下一个坐标的部分了。我该怎么做？

我想用来计算下一个坐标的牛顿法，是下面的方法：

function [zero,res,niter]=newton(f,df,x0,tol,nmax,varargin)
%NEWTON Find function zeros.
%   ZERO=NEWTON(FUN,DFUN,X0,TOL,NMAX) tries to find the zero ZERO of the 
%   continuous and differentiable function FUN nearest to X0 using the    Newton 
%   method. FUN and its derivative DFUN accept real scalar input x and returns 
%   a real scalar value. If the search fails an error message is displayed.
%   FUN and DFUN can also be inline objects.
%
%   [ZERO,RES,NITER]= NEWTON(FUN,...) returns the value of the residual in ZERO
%   and the iteration number at which ZERO was computed.

x = x0; 
fx = feval(f,x,varargin{:}); 
dfx = feval(df,x,varargin{:});
niter = 0; 
diff = tol+1;
resultVector = [x];
while diff >= tol && niter <= nmax
disp(x)
niter = niter + 1;
diff = - dfx\fx;
x = x + diff;     
diff = abs(diff);
fx = feval(f,x,varargin{:});        
dfx = feval(df,x,varargin{:});
resultVector = [resultVector x];
end
if niter > nmax
fprintf(['newton stopped without converging to the desired tolerance',...
'because the maximum number of iterations was reached\n']);
end
res = fx
zero = x;
disp('Number of iterations is:');
disp(niter);

plot(resultVector);
return

我想我没有完全理解'projecting all the original coordinates to coordinates in this subspace you can perform an m-1-dimensional Newton-Raphson procedure'部分。这是什么意思？我给牛顿法输入什么？

非常感谢！

Answer 1

我将回答如何在 Matlab 中通过采用牛顿法求解具有 n 个未知数的 n-1 方程组的问题。我的改编不是你通过研究发现的那种——它更简单。

牛顿法的思想是围绕某个猜测点将系统线性化，然后求解得到的线性系统。解决方案成为我们的下一个猜测。

如果你的系统是欠定的（方程比未知数少），线性系统也将是欠定的。 Matlab 的 \ 运算符可以很好地处理此类系统：returns 最小范数的解。那么，如果它有效，为什么不直接使用它呢？与具有 n 个未知数的 n 个方程的情况相比，不需要更改代码。

下面是我的实现：调用newtonund(1,1,1) 20步求根，即-1.226458955, 1.357877375, -0.031788652。我稍微更改了您的方程组，将 sin(z) 替换为 cos(z)。它也以原始形式工作，但结果（可以预见）是 0, 0, 0，这并不那么有趣。

function newtonund(x,y,z)
tol = 1e-6;
diff = tol+1;
n = 0
nmax = 1000;
disp('  n      x(n)            y(n)             z(n)          |f(x)|');
X = [x;y;z];
Y = F(X);
while diff>=tol & n<=nmax
    changeX = -dF(X)\Y;
    X = X + changeX;
    Y = F(X);
    diff = norm(changeX);
    n = n+1;
    fprintf('%3d %15.9f %15.9f %15.9f %10.5g \n', n, X(1), X(2), X(3), norm(Y));
end
if n>nmax
    disp('Failed to converge');
else
    disp('Root found');
end
end

function Y=F(X)
x=X(1); y=X(2); z=X(3);
Y=[x^3+y^2+z^2 ; x^2-y^3+cos(z)];
end

function J=dF(X)
x=X(1); y=X(2); z=X(3);
J=[3*x^2, 2*y, 2*z; 2*x, -3*y^2, -sin(z)];
end

changeX = -dF(X)\Y行表达了牛顿法。通常 dF(X) 是一个方阵；在不确定的情况下，它是矩形的，但该算法仍然有效。

如何在 Matlab 中使用雅可比和牛顿法求解（非线性）方程组

How to solve system of (non-linear) equations using Jacobian and Newton's Method in Matlab

matlab

newtons-method

numerical-methods

derivative

equation-solving