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

Question

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

    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, она теряет свое свойство как функцию, и тогда я не могу использовать ее в качестве аргумента во втором вызове для получения каждого элемента. Как это можно улучшить? Как я могу сделать программу / функцию подпрограммой, что ее вход является матрицей, подобной этой, и ее вывод будет ее производной???

Спасибо

Answer 1

Код не делает то, что вы думаете, что он делает. W это просто массив реалов. Вы просто заполняете записи с оценкой в ​​точке 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