在 Fortran 中计算雅可比矩阵
在牛顿法中,为了求解非线性方程组,我们需要找到雅可比矩阵和雅可比矩阵逆的行列式。
这是我的组件函数,
real function f1(x,y)
parameter (pi = 3.141592653589793)
f1 = log(abs(x-y**2)) - sin(x*y) - sin(pi)
end function f1
real function f2(x,y)
f2 = exp(x*y) + cos(x-y) - 2
end function f2
对于 2x2 的情况,我正在计算雅可比矩阵和雅可比矩阵逆的行列式,如下所示,
x = [2,2]
h = 0.00001
.
.
! calculate approximate partial derivative
! you can make it more accurate by reducing the
! value of h
j11 = (f1(x(1)+h,x(2))-f1(x(1),x(2)))/h
j12 = (f1(x(1),x(2)+h)-f1(x(1),x(2)))/h
j21 = (f2(x(1)+h,x(2))-f2(x(1),x(2)))/h
j22 = (f2(x(1),x(2)+h)-f2(x(1),x(2)))/h
! calculate the Jacobian
J(1,:) = [j11,j12]
J(2,:) = [j21,j22]
! calculate inverse Jacobian
inv_J(1,:) = [J(2,2),-J(1,2)]
inv_J(2,:) = [-J(2,1),J(1,1)]
DET=J(1,1)*J(2,2) - J(1,2)*J(2,1)
inv_J = inv_J/DET
.
.
我如何在 Fortran 中扩展它来评估在 n 点评估的 m 个函数的雅可比行列式?
In Newton's method, to solve a nonlinear system of equations we need to find the Jacobian matrix and the determinant of the inverse of the Jacobian matrix.
Here are my component functions,
real function f1(x,y)
parameter (pi = 3.141592653589793)
f1 = log(abs(x-y**2)) - sin(x*y) - sin(pi)
end function f1
real function f2(x,y)
f2 = exp(x*y) + cos(x-y) - 2
end function f2
For the 2x2 case I am computing the Jacobian matrix and determinant of the inverse of Jacobian matrix like this,
x = [2,2]
h = 0.00001
.
.
! calculate approximate partial derivative
! you can make it more accurate by reducing the
! value of h
j11 = (f1(x(1)+h,x(2))-f1(x(1),x(2)))/h
j12 = (f1(x(1),x(2)+h)-f1(x(1),x(2)))/h
j21 = (f2(x(1)+h,x(2))-f2(x(1),x(2)))/h
j22 = (f2(x(1),x(2)+h)-f2(x(1),x(2)))/h
! calculate the Jacobian
J(1,:) = [j11,j12]
J(2,:) = [j21,j22]
! calculate inverse Jacobian
inv_J(1,:) = [J(2,2),-J(1,2)]
inv_J(2,:) = [-J(2,1),J(1,1)]
DET=J(1,1)*J(2,2) - J(1,2)*J(2,1)
inv_J = inv_J/DET
.
.
How do I in Fortran extend this to evaluate a Jacobian for m functions evaluated at n points?
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

发布评论
评论(3)
关于评价我会分几个方面来回答。这很简单。您只需要一个点数组,如果这些点位于某个规则网格中,您甚至可能不需要它。
您可能有一个 x 数组和 y 数组,或者您可以有一个包含 x 和 y 分量的派生数据类型数组。
对于前者:
real, allocatable :: x(:), y(:)
x = [... !probably read from some data file
y = [...
do i = 1, size(x)
J(i) = Jacobian(f, x(i), y(i))
end do
如果你想同时拥有多个函数,问题总是出在函数的表示上。即使您有一个函数指针数组,您也需要手动对其进行编码。另一种方法是拥有一个完整的代数模块,您可以在其中输入一些表示函数的字符串,然后可以计算该函数,甚至可以符号地计算导数。这需要一个解析器、一个评估器,这是一项艰巨的任务。有一些库可以做到这一点。除非采取进一步的优化步骤(编译为机器代码),否则对此类导数的评估将会很慢。
导数的数值评估当然是可能的。它会稍微减慢收敛速度,具体取决于导数的近似阶数。您对数值导数进行两点之差。您可以根据多个点的值创建插值多项式以获得高阶近似(有限差分近似),但这会消耗机器周期。
通常我们可以使用@John Alexiou提到的自动差异工具。但在实践中,我更喜欢使用 MATLAB 分析求解雅可比行列式,然后使用其内置函数 fortran() 将结果转换为 f90 文件。以你的函数为例。只需将这些输入到 MATLAB 中
syms x y
Fval=sym(zeros(2,1));
Fval(1)=log(abs(x-y^2)) - sin(x*y) - sin(pi);
Fval(2)=exp(x*y) + cos(x-y) - 2;
X=[x;y];
Fjac=jacobian(Fval,X);
fortran(Fjac)
即可生成结果
Fjac(1,1) = -y*cos(x*y)-((-(x-y**2)/abs(-x+y**2)))/abs(-x+y**2)
Fjac(1,2) = -x*cos(x*y)+(y*((-(x-y**2)/abs(-x+y**2)))*2.0D0)/abs(-
&x+y**2)
Fjac(2,1) = -sin(x-y)+y*exp(x*y)
Fjac(2,2) = sin(x-y)+x*exp(x*y)
。您只需得到一个解析型雅可比 fortran 函数。
同时,由于秩不匹配,无法求解 mxn 矩阵的逆。您应该简化方程组以获得 nxn
雅各宾。
此外,当我们使用 Newton-Raphson 方法时,我们不求解 Jacobin 的逆,这对于大型系统来说既耗时又不准确。一个简单的方法是在 LAPACK 中使用 dgesv 来实现密集 Jacobin。因为我们只需要从线性方程组中求解向量x
Jx=-F
dgesv
使用LU分解和高斯消去法求解上述方程组,比求解逆矩阵要快得多。
如果方程组很大,可以使用UMFPACK
及其fortran接口模块mUMFPACK
来求解方程组,其中J是稀疏矩阵。或者在广泛使用的稀疏矩阵库 SPARSEKIT2
中使用子例程 ILUD
和 LUSOL
。
除此之外,还有大量其他方法试图更快、更准确地求解 Jx=-F,例如广义最小残差 (GMRES) 和稳定双-共轭梯度(BICGSTAB)是一篇文献。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
这是一个更灵活的雅可比计算器。
2×2 测试用例的结果正是您所期望的
的给定输入的符号计算检查结果
我根据
Console.f90
Here is a more flexible jacobian calculator.
The results with the 2×2 test case are what you expect
I check the results against a symbolic calculation for the given inputs of
Console.f90