计算科学 - 我如何计算这个函数的雅可比？ - 吾爱随笔录

我如何计算这个函数的雅可比？

计算科学雅可比

2021-12-16 04:43:48

我有两个向量 $r$ 和 $m$ . 两个向量都是 $N\times1$ .

一个函数计算为 -

$F(1:N) = \phi r + (r^3 + rm^2)$

$F(N+1:2N) = \phi m + (m^3+mr^2)$

F1 = (phi*r + r.^3 + r.*m.^2);
F2 = (phi*m + m.^3 + m.*r.^2);

F = [F1;F2];

我在计算这个函数的 Jacboian 时遇到了麻烦。这就是我所做的 -

J11 = phi*eye(N) + diag(3*r.^2 + m.^2);
J12 = diag(2*r.*m);
J21 = diag(2*m.*r);
J22 = phi*eye(N) + diag(3*m.^2 + r.^2);
J = [J11 J12; J21 J22];

其中 J 是雅可比行列式。我相信问题出在函数的一部分 $r$ 和 $m$ 彼此相乘，但我对此不是 100% 确定，即使我是，我也不确定如何解决这个问题。

上面的代码是用 MATLAB 编写的。任何帮助将非常感激。

1个回答

如果您使用了索引符号，您可能自己也注意到了这个错误，而我也很难看到它。以线路J11 = phi*eye(N) + (3*r.^2 + m.^2); 为例。您正在将矩阵phi*eye(N)和向量添加(3*r.^2 + m.^2)在一起。几年来，这是 MATLAB 中允许的操作，因为它在数据科学中有很好的用途，但是，我更喜欢以前抛出错误的行为，因为它会导致像这样的语义错误。例如，跑步[1 0; 0 1] + [2;2]返回

这对于数据科学来说可能是一种方便的操作，但对你来说却是一个错误的来源。

现在，让我通过用索引符号重写你的问题来指出这个问题。让 $F = [f_0, f_1, \dots, f_n]^T$ 在哪里 $f_i = \phi r_i + r_i^3 + r_i m_i^2$ 和 $G = [g_0, g_1, \dots, g_n]^T$ 在哪里 $g_i = \phi m_i + m_i^3 + m_i r_i^2$ . 我们要计算 $\nabla \begin{bmatrix} F \\ G \end{bmatrix}$ 这将是一个矩阵。让我们从写第一行开始，它由偏导数组成 $f_0 = \phi r_0 + r_0^3 + r_0 m_0^2$ ：

\partial_{r_{0}} f_{0} = ϕ + 3 r_{0}^{2} + m_{0}^{2},

$\partial_{r_0} f_0 = \phi + 3r_0^2 + m_0^2,$

\partial_{r_{i}} f_{0} = 0 for all i > 0,

$\partial_{r_i} f_0 = 0 \quad \text{for all } i>0,$

\partial_{m_{0}} f_{0} = 2 r_{0} m_{0},

$\partial_{m_0} f_0 = 2r_0m_0,$ 最后

\partial_{m_{i}} f_{0} = 0 for all i > 0.

$\partial_{m_i} f_0 = 0 \quad \text{for all } i>0.$ 因此，鉴于我正确解释了您的符号，所有雅可比计算都是错误的。（但我不知道如何理解

F

$F$ 否则）

因此，雅可比应该计算为

J11 = phi*eye(N) + diag(3*r.^2 + m.^2);
J12 = diag(2*r.*m);
J21 = diag(2*m.*r);
J22 = phi*eye(N) + diag(3*m.^2 + r.^2);

J = [J11 J12; J21 J22];

根据未知数的数量，您可能希望使用 MATLAB 提供的稀疏数据结构（speye、spdiags 等）。

其它你可能感兴趣的问题

上一篇gmres MATLAB中的relres 下一篇如何从 2D FFT 幅度与波数图获得准确的波长值，就像从 1D FFT 幅度与波数图获得的一样？