计算科学 - 求解这个非线性方程组 - 吾爱随笔录

求解这个非线性方程组

计算科学牛顿法非线性方程

2021-12-01 05:14:02

假设我有这组方程：

a = x + z (1)

$a = x + z\qquad (1)$

b = y + \frac{z}{2} (2)

$b = y + \frac{z}{2}\qquad (2)$

z = k_{0} x \sqrt{y} (3)

$z = k_0x\sqrt{y}\qquad (3)$

其中，和。和的值假设 x 和 y 都不能小于零，我如何稳健地解决这个系统？我尝试了 Newton-Raphson 方法并得到了以下残差和雅可比： $a$ $b$ $\in \mathbb{R}$ $k_0 > 0$ $a$ $b$

R (x, y) = {\begin{matrix} x + z - a \\ y + \frac{z}{2} - b \\ k_{0} x \sqrt{y} - z \end{matrix}}

$R(x,y) = \begin{Bmatrix}x+z-a \\y+\frac{z}{2}-b \\ k_0x\sqrt{y} - z\end{Bmatrix}$

J (x, y) = [\begin{matrix} 1 & 0 & 1 \\ 0 & 1 & \frac{1}{2} \\ \frac{k_{0} \sqrt{y}}{2} & \frac{k_{0} x}{2 \sqrt{y}} & - 1 \end{matrix}]

$J(x,y) = \begin{bmatrix}1 && 0 && 1 \\ 0 && 1 && \frac{1}{2} \\ \frac{k_0\sqrt{y}}{2} && \frac{k_0x}{2\sqrt{y}} && -1 \end{bmatrix}$

这是MATLAB代码：

% User-defined inputs
a = 1;
b = 0;
k0 = 10;

% Initial guess
%u = [1;1;1];
u = rand(3,1)*2;

% Newton-Raphson solver
i = 0; max_iters = 30; TOL = 1e-12;
fprintf('a = %f\nb = %f\nk0 = %f\nInitial guess: (%f,%f,%f)\n\n\tResidual_x\tResidual_y\tResidual_z\tResidual_rel\n',a,b,k0,u(1),u(2),u(3));
while (true)

    % Residual function
    F = [u(1)+u(3)-a;u(2)+0.5*u(3)-b;k0*u(1)*sqrt(u(2))-u(3)];

    % Initial residual
    if i == 0
        R_init = norm(F);
    % Maximum iterations
    elseif i == max_iters
        fprintf('\tSolution did not converge.\n');
        break;
    end

    % Absolute residual
    R = norm(F);

    % Relative residual
    R_rel = R/R_init;
    fprintf('\t%e\t%e\t%e\t%e)\n',norm(F(1)),norm(F(2)),norm(F(3)),R_rel);

    % Check for convergence
    if (R_rel < TOL)
        fprintf('\tConverged.\n');
        break;
    end

    % Jacobian matrix
    J = [1,0,1;0,1,0.5;u(1)*k0*sqrt(u(2)),0.5*k0*u(1)/sqrt(u(2)),-1];

    % Solve and update
    u = u-J\F;
    i = i + 1;
end
fprintf('\nF

似乎这并不总是有效。例如，如果在 iterate我有，则求解器会爆炸（因为会导致错误）。最初的猜测似乎也发挥了重要作用，因为如果您使用随机猜测运行上述代码，每次都会得到不同的解决方案。 $i$ $y^i = 0$ $\frac{1}{\sqrt{0}}$

也就是说，我还能如何（或更好地）解决我上面的问题？或者如果我坚持这个 Newton-Raphson 方案，我应该如何选择我的初始猜测？

1个回答

在 sage 或任何其他 CAS 的帮助下，这可以很容易地简化为一系列单变量方程。（我替换以简化事情 $\sqrt{y}\mapsto y$ ，唯一改变的是下面只有的解决方案才有效。） $y\geq0$

R = PolynomialRing(QQ, 'a,b,k,x,y,z', order='invlex')
a,b,k,x,y,z = R.gens()
I = R.ideal(x+z-a, y*y+z/2-b, k*x*y-z)
# should be 3, meaning finitely many solutions for each set of parameters
print(I.dimension()) 
I.groebner_basis()

\begin{aligned} [ & z + x - a, \\ y^{2} - \frac{1}{2} x - b + \frac{1}{2} a, \\ x y - a y + \frac{1}{2} k x^{2} + b k x - \frac{1}{2} a k x, \\ a k y - \frac{1}{2} k^{2} x^{2} - b k^{2} x + \frac{1}{2} a k^{2} x + x - a, \\ k^{2} x^{3} + 2 b k^{2} x^{2} - a k^{2} x^{2} - 2 x^{2} + 4 a x - 2 a^{2}] \end{aligned}

$\begin{aligned} \big[&z + x - a, \\&y^{2} - \tfrac{1}{2} x - b + \tfrac{1}{2} a, \\&x y - a y + \tfrac{1}{2} k x^{2} + b k x - \frac{1}{2} a k x, \\&a k y - \frac{1}{2}k^{2} x^{2} - b k^{2} x + \tfrac{1}{2} a k^{2} x + x - a, \\&k^{2} x^{3} + 2 b k^{2} x^{2} - a k^{2} x^{2} - 2 x^{2} + 4 a x - 2 a^{2}\big] \end{aligned}$ 因此必须满足三次方程（另一个要尝试的命令是。）

x

$x$

- 2 a^{2} + 4 a x - (2 + a k^{2} - 2 b k^{2}) x^{2} + k^{2} x^{3} = 0

$-2 a^2+4 a x-(2+a k^2-2 b k^2) x^2+k^2 x^3 = 0$ I.elimination_ideal([y,z])

根据最多有三个可能的值。对于的每个可能值，您可以求解，然后求解线性方程，最多给出三个不同的解。然后可以检查每个解决方案的有效性并在必要时丢弃。 $a,k$ $x$ $x$ $a=x+z$ $y=z/(kx)$

其它你可能感兴趣的问题

上一篇如何在 MATLAB 中定义微分函数中的代数方程？下一篇Lighthill 源项的有限差分