Я пытаюсь численно выполнить 2D-интеграл, используя библиотеку quadpack double . Поскольку на этом веб-сайте нет встроенного движка Latex, я прилагаю картинку с математикой.
Также ниже прикрепляю свою реализацию этой операции в модуле test_module.
module test_module
use f90_kind ! defines dp
use physconst ! defines pi
use global_variables ! defines the variables for quadpack
use csv_file ! allows me to write to a csv file
implicit none
real(dp) :: a, b, c ! parameters I want to set
contains
real(dp) function multi_val_func(theta)
real(dp), intent(in) :: theta
multi_val_func = exp(theta**2 + c*theta)
end function multi_val_func
real(dp) function theta_integral(phi)
real(dp), intent(in) :: phi
c = phi*(a+b)
! Quadpack variables
epsabs = 1.0E-14_dp
epsrel = 0.0E0_dp
key = 6
ilow = 0.0E0_dp
ihigh = 2.0E0_dp*pi
call dqage(multi_val_func, ilow, ihigh, epsabs, epsrel, key, limit, &
res, abserr, neval, ier, alist, blist, rlist, elist, iord, last)
theta_integral = res
end function theta_integral
subroutine phi_integral(final_result)
real(dp),intent(out) :: final_result
a = 1.0E0_dp
b = 1.0E0_dp
! Quadpack variables
epsabs = 1.0E-14_dp
epsrel = 0.0E0_dp
key = 6
ilow = 0.0E0_dp
ihigh = pi
call dqage(theta_integral, ilow, ihigh, epsabs, epsrel, key, limit, &
res, abserr, neval, ier, alist, blist, rlist, elist, iord, last)
final_result = res
print *, final_result
end subroutine phi_integral
end module test_module
В основном файле я определяю переменную result, а затем call phi_integral(result). Я компилирую этот код, используя следующий make-файл
FFLAGS = -O0 -fcheck=all -ffree-line-length-none
debug:
gfortran -c $(FFLAGS) $(srcdir)/f90_kind.f90
gfortran -c $(FFLAGS) $(srcdir)/physconst.f90
gfortran -c $(FFLAGS) $(srcdir)/global_variables.f90
gfortran -c $(FFLAGS) $(extlib_qp)/quadpack_double.f90
gfortran -c $(FFLAGS) $(extlib_csv)/csv_file.f90
gfortran -c $(FFLAGS) $(srcdir)/multi_var_integration.f90
gfortran -c $(FFLAGS) $(srcdir)/main.f90
gfortran *.o -o debug -lblas -llapack
rm -f *.o *.mod
Код компилируется, но когда я его запускаю, я получаю следующую ошибку:
Fortran runtime error: Recursive call to nonrecursive procedure 'dqage'
Error termination. Backtrace:
#0 0x103a28d3d
#1 0x103a299f5
#2 0x103a29fd6
#3 0x1039f94b6
#4 0x1039dee59
#5 0x1039fad68
#6 0x1039f9948
#7 0x1039dec7d
#8 0x1039def61
#9 0x1039defa3
Любое предложение по решению этой проблемы приветствуется.
Не каждая подпрограмма или функция в Фортране является реентерабельной или, как ее называет Фортран, «рекурсивной». До Fortran 2018 подпрограмма, которую можно вызвать, находясь уже внутри себя (пусть даже косвенно), должна быть отмечена атрибутом recursive.
Даже если подпрограмма помечена таким образом, она может работать не так, как ожидалось, если она содержит статические (сохраненные) данные или если она использует глобальные переменные (общие блоки, модули).
Я не могу сразу увидеть какие-либо такие статические данные, используемые в http://www.netlib.no/netlib/quadpack/dqage.f, поэтому они могут работать после добавления recursive.
Однако, как указал Франческалус, вам лучше иметь специальную библиотеку 2d-квадратур.
Интеграл theta можно выполнить с помощью функции ошибок erf. См. вольфрам альфа.
Результатом этого является функция одной переменной phi, которую можно вычислить через quadpack. Таким образом, вам не нужна многомерная интегральная процедура.
Вы понимаете, что означает «Рекурсивный вызов нерекурсивной процедуры» и каковы последствия? Вероятно, вам лучше попытаться использовать квадратурную процедуру/пакет, предназначенный для двумерных интегралов.