Как получить матрицу, функции которой являются ее элементами?

Question

Я использую fortran, и я пытаюсь взять производную от матрицы, которая является ее элементами.

program derivada_matrix 

    integer, parameter :: matrix_size = 5 
    integer :: i,j 
    real(8) :: time = 1.0 
    real(8),dimension (matrix_size, matrix_size) :: W 
    real(8),dimension (matrix_size, matrix_size) :: dW 

    call potent(time,W) 
    do i = 1, matrix_size 
     do j=1, matrix_size 
     call Derivada(time,W(i,j),dW(i,j)) 
     end do 
    end do 

    print*, 'matrix' 
    print*, W 
    print*, 'derivada', dW 


    end program 


    Subroutine Derivada (x1,F,D) 
    implicit none 
    Real*8 :: x1 
    Real*8 :: h= 1.0E-6 
    integer, parameter :: matrix_size = 5 
    real*8 :: D,F 
    external F 
    D = (1.0*F(x1-2*h) - 8.0*F(x1-h) + 8.0*F(x1+h) - 1.0*F(x1+2*h))/(12.0*h) 

return 
End subroutine Derivada 

subroutine potent(T,W) 
implicit none 
integer, parameter :: matrix_size = 5 
real(8),dimension(matrix_size,matrix_size) :: W 
Real(8):: T 
integer :: i,j 
do i = 1, matrix_size 
    do j=i,matrix_size 
     W(i,j) = 0.0 
     W(j,i) = W(i,j) 
    end do 
    W(i,i) = cos(T) 
end do 
RETURN 
END subroutine potent 

В основном первая подпрограмма создает тестовую матрицу с функцией (косинуса) на главной диагонали и нули в других местах, а второй подпрограммы он должен его производного. Это сообщение об ошибке/предупреждение я получаю

 call Derivada(time,W(i,j),dW(i,j)) 
         1 
    Warning: Expected a procedure for argument ‘f’ at (1) 

Ошибка/предупреждение я получаю во втором вызове. Думаю, потому что, когда я создаю матрицу W, она утрачивает свое свойство как функцию, а затем я не могу использовать в качестве аргумента во втором вызове для вывода каждого элемента. Как это можно улучшить? Как я могу сделать подпрограмму программы/функции, чтобы ее вход был такой матрицей, и ее выход был бы ее производной ???

Thanks

Answer 1

Код не подходит к тому, что вы думаете. W - всего лишь массив реалов. Вы просто заполняете записи cos, оцениваемыми в точке T. Это не набор функций.

Вот более современный подход к тому, что вы хотите сделать. Вместо использования extern-функций он использует указатели на процедуры и обертывает их производным типом, поэтому вы можете создать массив из них. Затем вы указываете на функции, которые вы хотите, чтобы взять производные.

Ограничение этого подхода заключается в том, что вы должны соблюдать интерфейсы, в этом случае он будет работать только с функциями, которые принимают один двойной аргумент. Вы можете развернуть это несколькими способами, например, используя необязательные аргументы, обернуть их производным типом или использовать несколько интерфейсов, но это зависит от вас.

module deriv 
    use, intrinsic :: iso_fortran_env, only: dp => real64 
    implicit none 

    ! Wrapper for a procedure pointer so we can make an array of them 
    type ptr_wrapper 
    procedure(f), nopass, pointer :: func 
    end type ptr_wrapper 


    ! Use an interface for the function rather than "extern" functions 
    ! Change the interface to whatever you want. You could take multiple reals, integers, etc 
    abstract interface 
    pure function f(x1) 
     import 
     real(dp), intent(in) :: x1 
     real(dp) :: f 
    end function f 
    end interface 


contains 

    function derivada(x1, fx) result(d) 
    implicit none 
    real(dp), intent(in) :: x1 
    procedure(f), pointer :: fx 
    real(dp) :: d 

    real(dp) :: h = 1.0E-6 

    d = (1.0*fx(x1-2*h) - 8.0*fx(x1-h) + 8.0*fx(x1+h) - 1.0*fx(x1+2*h))/(12.0*h) 

    end function derivada 


    pure function test_func(x1) 
    real(dp), intent(in) :: x1 
    real(dp) test_func 

    test_func = 2d0 * x1 
    end function test_func 

    pure function test_func2(x1) 
    real(dp), intent(in) :: x1 
    real(dp) test_func2 

    test_func2 = x1**2 
    end function test_func2 


end module deriv 

program derivada_matrix 
    use, intrinsic :: iso_fortran_env, only: dp => real64 
    use deriv 
    implicit none 

    integer i, j 
    real(dp) :: dx 


    type(ptr_wrapper) :: w(2,2) 

    !Point to the function that you want each index to point to 
    w(1,1)%func => test_func 
    w(1,2)%func => test_func2 
    ! Pointing to the generic function wasn't working so I had to point to the specific double value trig functions 
    w(2,1)%func => dcos 
    w(2,2)%func => dsin 

    do i = 1, 2 
    do j = 1, 2 
     dx = derivada(1d0, w(i,j)%func) 
     print *, dx 
    end do 
    end do 


end program derivada_matrix 

И выводит производную этих функций при х = 1,0

2.00000000000000  
    2.00000000011102  
-0.841470984776962  
    0.540302305853706