在 Fortran 中计算雅可比矩阵

发布于 01-15 14:30 字数 1062 浏览 4 评论 0原文

在牛顿法中,为了求解非线性方程组,我们需要找到雅可比矩阵和雅可比矩阵逆的行列式。

这是我的组件函数,

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 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(3

嗫嚅2025-01-22 14:30:23

这是一个更灵活的雅可比计算器。

2×2 测试用例的结果正是您所期望的

arguments (x)
   2.00000000000000
   2.00000000000000

 values (y)
   1.44994967586787
   53.5981500331442

 Jacobian
  0.807287239448229        3.30728724371454
   109.196300248300        109.196300248300

的给定输入的符号计算检查结果

我根据fig2

Console.f90

program Console1
use ISO_FORTRAN_ENV
implicit none

! Variables
integer, parameter :: wp = real64
real(wp), parameter :: pi = 3.141592653589793d0

! Interfaces
interface 
    function fun(x,n,m) result(y)
    import
        integer, intent(in) :: n,m
        real(wp), intent(in) :: x(m)
        real(wp) :: y(n)
    end function
end interface

real(wp) :: h
real(wp), allocatable :: x(:), y(:), J(:,:)
! Body of Console1

x = [2d0, 2d0]
h = 0.0001d0

print *, "arguments"
print *, x(1)
print *, x(2)

y = test(x,2,2)
print *, "values"
print *, y(1)
print *, y(2)

J = jacobian(test,x,2,h)

print *, "Jacobian"
print *, J(1,:)
print *, J(2,:)

contains

function test(x,n,m) result(y)
! Test case per original question
    integer, intent(in) :: n,m
    real(wp), intent(in) :: x(m)
    real(wp) :: y(n)
    
    y(1) = log(abs(x(1)-x(2)**2)) - sin(x(1)*x(2)) - sin(pi)  
    y(2) = exp(x(1)*x(2)) + cos(x(1)-x(2)) - 2 
    
end function
   
function jacobian(f,x,n,h) result(u)
    procedure(fun), pointer, intent(in) :: f
    real(wp), allocatable, intent(in) :: x(:)
    integer, intent(in) :: n
    real(wp), intent(in) :: h
    real(wp), allocatable :: u(:,:)
    integer :: j, m
    real(wp), allocatable :: y1(:), y2(:), e(:)
    
    m = size(x)
    allocate(u(n,m))
    
    do j=1, m
        e = element(j, m)    ! Get kronecker delta for j-th value
        y1 = f(x-e*h/2,n,m)
        y2 = f(x+e*h/2,n,m)
        u(:,j) = (y2-y1)/h   ! Finite difference for each column
    end do        
    
end function

function element(i,n) result(e)
! Kronecker delta vector. All zeros, except the i-th value.
integer, intent(in) :: i, n
real(wp) :: e(n)
    e(:) = 0d0
    e(i) = 1d0        
end function

end program Console1

Here is a more flexible jacobian calculator.

The results with the 2×2 test case are what you expect

arguments (x)
   2.00000000000000
   2.00000000000000

 values (y)
   1.44994967586787
   53.5981500331442

 Jacobian
  0.807287239448229        3.30728724371454
   109.196300248300        109.196300248300

I check the results against a symbolic calculation for the given inputs of

fig2

Console.f90

program Console1
use ISO_FORTRAN_ENV
implicit none

! Variables
integer, parameter :: wp = real64
real(wp), parameter :: pi = 3.141592653589793d0

! Interfaces
interface 
    function fun(x,n,m) result(y)
    import
        integer, intent(in) :: n,m
        real(wp), intent(in) :: x(m)
        real(wp) :: y(n)
    end function
end interface

real(wp) :: h
real(wp), allocatable :: x(:), y(:), J(:,:)
! Body of Console1

x = [2d0, 2d0]
h = 0.0001d0

print *, "arguments"
print *, x(1)
print *, x(2)

y = test(x,2,2)
print *, "values"
print *, y(1)
print *, y(2)

J = jacobian(test,x,2,h)

print *, "Jacobian"
print *, J(1,:)
print *, J(2,:)

contains

function test(x,n,m) result(y)
! Test case per original question
    integer, intent(in) :: n,m
    real(wp), intent(in) :: x(m)
    real(wp) :: y(n)
    
    y(1) = log(abs(x(1)-x(2)**2)) - sin(x(1)*x(2)) - sin(pi)  
    y(2) = exp(x(1)*x(2)) + cos(x(1)-x(2)) - 2 
    
end function
   
function jacobian(f,x,n,h) result(u)
    procedure(fun), pointer, intent(in) :: f
    real(wp), allocatable, intent(in) :: x(:)
    integer, intent(in) :: n
    real(wp), intent(in) :: h
    real(wp), allocatable :: u(:,:)
    integer :: j, m
    real(wp), allocatable :: y1(:), y2(:), e(:)
    
    m = size(x)
    allocate(u(n,m))
    
    do j=1, m
        e = element(j, m)    ! Get kronecker delta for j-th value
        y1 = f(x-e*h/2,n,m)
        y2 = f(x+e*h/2,n,m)
        u(:,j) = (y2-y1)/h   ! Finite difference for each column
    end do        
    
end function

function element(i,n) result(e)
! Kronecker delta vector. All zeros, except the i-th value.
integer, intent(in) :: i, n
real(wp) :: e(n)
    e(:) = 0d0
    e(i) = 1d0        
end function

end program Console1
-黛色若梦2025-01-22 14:30:23

关于评价我会分几个方面来回答。这很简单。您只需要一个点数组,如果这些点位于某个规则网格中,您甚至可能不需要它。

您可能有一个 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

如果你想同时拥有多个函数,问题总是出在函数的表示上。即使您有一个函数指针数组,您也需要手动对其进行编码。另一种方法是拥有一个完整的代数模块,您可以在其中输入一些表示函数的字符串,然后可以计算该函数,甚至可以符号地计算导数。这需要一个解析器、一个评估器,这是一项艰巨的任务。有一些库可以做到这一点。除非采取进一步的优化步骤(编译为机器代码),否则对此类导数的评估将会很慢。

导数的数值评估当然是可能的。它会稍微减慢收敛速度,具体取决于导数的近似阶数。您对数值导数进行两点之差。您可以根据多个点的值创建插值多项式以获得高阶近似(有限差分近似),但这会消耗机器周期。

I will answer about evaluation in different points. This is quite simple. You just need an array of points, and if the points are in some regular grid, you may not even need that.

You may have an array of xs and array of ys or you can have an array of derived datatype with x and y components.

For the former:

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

If you want to have many functions at the same time, the problem is always in representing functions. Even if you have an array of function pointers, you need to code them manually. A different approach is to have a full algebra module, where you enter some string representing a function and you can evaluate such function and even compute derivatives symbolically. That requires a parser, an evaluator, it is a large task. There are libraries for this. Evaluation of such a derivative will be slow unless further optimizing steps (compiling to machine code) are undertaken.

Numerical evaluation of the derivative is certainly possible. It will slow the convergence somewhat, depending on the order of the approximation of the derivative. You do a difference of two points for the numerical derivative. You can make an interpolating polynomial from values in multiple points to get a higher-order approximation (finite difference approximations), but that costs machine cycles.

朮生2025-01-22 14:30:23

通常我们可以使用@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 中使用子例程 ILUDLUSOL

除此之外,还有大量其他方法试图更快、更准确地求解 Jx=-F,例如广义最小残差 (GMRES) 和稳定双-共轭梯度(BICGSTAB)是一篇文献。

Normally we can use auto difference tools as @John Alexiou mentioned. However in practise I prefer using MATLAB to analytically solve out the Jacobian and then use its build-in function fortran() to convert the result to a f90 file. Take your function as an example. Just type these into 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) 

which will yield

  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)

to you. You just get an analytical Jacobian fortran function.

Meanwhile, it is impossible to solve the inverse of a mxn matrix because of rank mismatching. You should simplify the system of equations to get a nxn Jacobin.

Additionally, when we use Newton-Raphson's method we do not solve the inverse of the Jacobin which is time-consuming and inaccurate for a large system. An easy way is to use dgesv in LAPACK for dense Jacobin. As we only need to solve the vector x from system of linear equations

Jx=-F

dgesv use LU decomposition and Gaussian elimination to solve above system of equations which is extremely faster than solving inverse matrix.

If the system of equations is large, you can use UMFPACK and its fortran interface module mUMFPACK to solve the system of equations in which J is a sparse matrix. Or use subroutine ILUD and LUSOL in a wide-spread sparse matrix library SPARSEKIT2.

In addition to these, there are tons of other methods which try to solve the Jx=-F faster and more accurate such as Generalized Minimal Residual (GMRES) and Stabilized Bi-Conjugate Gradient (BICGSTAB) which is a strand of literature.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文