计算科学 - 关于在寻根算法中确定雅可比的困惑 - 吾爱随笔录

关于在寻根算法中确定雅可比的困惑

计算科学 Python scipy 根雅可比寻根

2021-12-18 04:05:00

我编写了一些 Python 代码来确定以下非线性方程的数值根：

f_{m} = \tan λ_{m} - \frac{λ_{m}}{1 + a}

$f_m=\tan\lambda_m - \frac{\lambda_m}{1+a}$ 在哪里

λ_{m} > 0

$\lambda_m\gt0$ 和

a \geq 0

$a\geq0$ . 代码是：

import numpy as np
from scipy.optimize import root

func = lambda λm, a: np.tan(λm) - λm/(1+a)
λm0 = np.linspace(np.pi, m*np.pi, m) # initial guess
sol = root(func, λm0, args=(a), method='lm') # least-square optimization

这很好用，例如对于输入结果是：这是的根，可以很容易地检查。 $m, a = 3, 1$

λ_{m} = [4.27478227, 7.59654602, 10.81267333]

$\lambda_m=[4.27478227, 7.59654602, 10.81267333]$

f_{m} = 0

$f_m=0$

m = 1, 2, 3

$m=1,2,3$

作为下一步，我想提供雅可比行列式，我认为应该是：即矩阵形式具有上述值：

J_{m, n} = \frac{\partial f_{m}}{\partial λ_{n}} = {\begin{cases} \frac{1}{\cos^{2} λ_{m}} - \frac{1}{1 + a} & m = n \\ 0 & m \neq n \end{cases}

$J_{m,n}=\frac{\partial f_{m}}{\partial\lambda_{n}}=\begin{cases} \frac{1}{\cos^2\lambda_{m}}-\frac{1}{1+a} & m=n\\ 0 & m\neq n \end{cases}$

λ_{m}

$\lambda_m$

J = [\begin{matrix} 5.06844087 & 0 & 0 \\ 0 & 14.92687786 & 0 \\ 0 & 0 & 29.72847616 \end{matrix}]

$\boldsymbol{J}=\begin{bmatrix} 5.06844087 & 0 & 0\\ 0 & 14.92687786 & 0\\ 0 & 0 & 29.72847616 \end{bmatrix}$

然而，优化例程输出的近似雅可比矩阵是：

\tilde{J} = [\begin{matrix} - 29.72850523 & 0 & 1 \\ 0 & - 14.92688569 & 0 \\ 0 & 0 & 5.06844398 \end{matrix}]

$\boldsymbol{\tilde{J}}=\begin{bmatrix} -29.72850523 & 0 & 1\\ 0 & -14.92688569 & 0\\ 0 & 0 & 5.06844398 \end{bmatrix}$

显然，对角线上的近似导数的大小对应于我对雅可比行列式的分析值，但这些值的顺序有些不同（从向后运行），并且一些值是负数。此外，存在一个额外的（这些虚假的似乎随着的增长而出现更多），它不应该存在。 $n=3,2,1$ $\tilde{J}_{1,3}=1$ $n$

有人可以解释近似雅可比行列式和解析雅可比行列式之间的区别吗？

1个回答

如果您查看脚本的完整输出，sol则有字段ipvt和qtf. 这两个字段都来自MINPACK（http://www.netlib.org/minpack/lmdif.f），这是这里的后端（https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize /_root.py#L246 ; https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.leastsq.html）。

首先，ipvt解释了变量排序发生了什么——在处理雅可比行列时有一个旋转步骤，当我运行你的代码时，枢轴是，与矩阵的排列相匹配。第二个，在文档中给出，在我的例子中是. 这解释了前两个条目的符号变化，但幅度 ——我不知道为什么。 $[3,2,1]$ qft(q transpose)*fvec[ -2.52376384e-07, -1.46294147e-07, 5.03946448e-07] $10^7$

我注意到的另一件事是它lmdif不会产生实际的雅可比行列式 - 请注意该字段被称为fjac（如在 lmdif 中），而不是jac（如在 scipy 中）。引用：

c       fjac is an output m by n array. the upper n by n submatrix
c         of fjac contains an upper triangular matrix r with
c         diagonal elements of nonincreasing magnitude such that
c
c                t     t           t
c               p *(jac *jac)*p = r *r,
c
c         where p is a permutation matrix and jac is the final
c         calculated jacobian. column j of p is column ipvt(j)
c         (see below) of the identity matrix. the lower trapezoidal
c         part of fjac contains information generated during
c         the computation of r.

我没有对此进行测试，但是如果您尝试使用耦合方程组，可能会出现更多的不一致。

其它你可能感兴趣的问题

上一篇这个非线性对流-扩散方程可以离散化为只需要求解 SPD 系统吗？下一篇MATLAB 中的 3D 可视化