f2py, função Python que retorna uma matriz (função com valor de vetor)
No seguinte Python eu tenho cinco funções contidas na matriz retornada porfunc
que eu tenho que integrar. O código chama um módulo externo de Fortran gerado usandof2py
:
import numpy as np
from numpy import cos, sin , exp
from trapzdv import trapzdv
def func(x):
return np.array([x**2, x**3, cos(x), sin(x), exp(x)])
if __name__ == '__main__':
xs = np.linspace(0.,20.,100)
ans = trapzdv(func,xs,5)
print 'from Fortran:', ans
print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.), exp(20.)])
A rotina do Fortran é:
subroutine trapzdv(f,xs,nf,nxs,result)
integer :: I
double precision :: x1,x2
integer, intent(in) :: nf, nxs
double precision, dimension(nf) :: fx1,fx2
double precision, intent(in), dimension(nxs) :: xs
double precision, intent(out), dimension(nf) :: result
external :: f
result = 0.0
do I = 2,nxs
x1 = xs(I-1)
x2 = xs(I)
fx1 = f(x1)
fx2 = f(x2)
result = result + (fx1+fx2)*(x2-x1)/2
enddo
return
end
O problema é que Fortran está apenas integrando a primeira função emfunc(x)
. Veja o resultado da impressão:
from Fortran: [ 2666.80270721 2666.80270721 2666.80270721 2666.80270721 2666.80270721]
exact: [ 2.66666667e+03 4.00000000e+04 9.12945251e-01 -4.08082062e-01 4.85165195e+08]
Uma maneira de trabalhar isso é modificarfunc(x)
para retornar o valor de uma determinada posição na matriz de funções:
def func(x,i):
return np.array([x**2, x**3, cos(x), sin(x), exp(x)])[i-1]
E então mude a rotina Fortran para chamar a função com dois parâmetros:
subroutine trapzdv(f,xs,nf,nxs,result)
integer :: I
double precision :: x1,x2,fx1,fx2
integer, intent(in) :: nf, nxs
double precision, intent(in), dimension(nxs) :: xs
double precision, intent(out), dimension(nf) :: result
external :: f
result = 0.0
do I = 2,nxs
x1 = xs(I-1)
x2 = xs(I)
do J = 1,nf
fx1 = f(x1,J)
fx2 = f(x2,J)
result(J) = result(J) + (fx1+fx2)*(x2-x1)/2
enddo
enddo
return
end
O que funciona:
from Fortran: [ 2.66680271e+03 4.00040812e+04 9.09838195e-01 5.89903440e-01 4.86814128e+08]
exact: [ 2.66666667e+03 4.00000000e+04 9.12945251e-01 -4.08082062e-01 4.85165195e+08]
Mas aquifunc
é chamado 5 vezes mais do que o necessário (no caso realfunc
tem mais de 300 funções, por isso será chamado 300 vezes mais do que o necessário).
func(x)
? Em outras palavras, faça a construção de Fortranfx1 = f(x1)
como uma matriz com 5 elementos correspondentes às funções emfunc(x)
.OBS: Estou compilando usandof2py -c --compiler=mingw32 -m trapzdv trapzdv.f90