在下面的Python中,我有五个函数包含在func
返回的数组中,我必须集成这些函数。代码调用使用f2py
生成的外部Fortran模块:
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.)])
Fortran例程是:
^{pr2}$问题是Fortran只对func(x)
中的第一个函数进行积分。
查看打印结果:
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]
一种解决方法是修改func(x)
以返回给定的
在函数数组中的位置:
def func(x,i):
return np.array([x**2, x**3, cos(x), sin(x), exp(x)])[i-1]
然后更改Fortran例程,用两个参数调用函数:
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
哪个有效:
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]
但是这里func
的调用比必要的多5倍(在实际情况中func
有超过300个函数,因此它将被调用300倍以上)。在
func(x)
返回的所有数组?换句话说,将Fortran build fx1 = f(x1)
作为一个数组,其中5个元素对应于func(x)
中的函数。在OBS:我正在使用f2py -c --compiler=mingw32 -m trapzdv trapzdv.f90
进行编译
不幸的是,不能将数组从python函数返回到Fortran中。为此您需要一个子程序(这意味着它是用
call
语句调用的),而这是f2py
不允许的。在在Fortran 90中,您可以创建返回数组的函数,但这又不是
f2py
所能做到的,特别是因为您的函数不是Fortran函数。在唯一的选择是使用循环解决方案,或者重新设计python和Fortran的交互方式。在
尽管这个答案并不能解决这个问题,但在Cython中它也是一个解决方法。这里对向量值函数实现了梯形规则和多项式积分器。我将下面的代码放入
integratev.pyx
:采用以下试验:
^{pr2}$给予:
多边形可以是任何顺序,这可能会给大多数情况下更好的结果。。。在
相关问题 更多 >
编程相关推荐