5

在下面的 Python 中，我有五个函数包含在返回的数组中func，我必须集成这些函数。代码调用使用生成的外部 Fortran 模块f2py：

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 例程是：

      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

问题是 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次以上）。

  • 有谁知道更好的解决方案让 Fortran 识别返回的所有数组func(x)？换句话说，将 Fortran 构建fx1 = f(x1)为一个包含 5 个元素的数组，这些元素对应于func(x).

OBS：我正在编译使用f2py -c --compiler=mingw32 -m trapzdv trapzdv.f90

4

2 回答 2

2

不幸的是，您不能将数组从 python 函数返回到 Fortran。为此，您需要一个子例程（意味着使用call语句调用它），这是f2py不允许您做的事情。

在 Fortran 90 中，您可以创建返回数组的函数，但这又不是f2py可以做到的，特别是因为您的函数不是 Fortran 函数。

您唯一的选择是使用循环解决方法，或者重新设计您希望 python 和 Fortran 交互的方式。

于 2013-07-04T18:26:02.480 回答
1

尽管这个答案不能解决问题，但它是在 Cython 中做同样的解决方法。这里梯形规则和多项式积分器是针对向量值函数实现的。下面的代码我放在一个integratev.pyx：

import numpy as np
from numpy.linalg import inv
cimport numpy as np
FLOAT = np.float32
ctypedef np.float_t FLOAT_t

def trapzv(f, np.ndarray xs, int nf):
    cdef int nxs = xs.shape[0]
    cdef np.ndarray ans = np.zeros(nf, dtype=FLOAT)
    cdef double x1, x2
    for i in range(1,nxs):
        x1 = xs[i-1]
        x2 = xs[i]
        ans += (f(x2)+f(x1))*(x2-x1)/2.
    return ans

def poly(f, np.ndarray xs, int nf, int order=2):
    cdef int nxs = xs.shape[0]
    cdef np.ndarray ans = np.zeros(nf, dtype=FLOAT)
    cdef np.ndarray xis = np.zeros(order+1, dtype=FLOAT)
    cdef np.ndarray ais
    if nxs % (order+1) != 0:
        raise ValueError("poly: The size of xs must be a multiple of 'order+1'")
    for i in range(order,nxs,order):
        xis = xs[i-order:i+1]
        X = np.concatenate([(xis**i)[:,None] for i in range(order+1)], axis=1)
        ais = np.dot( inv(X), f(xis).transpose() )
        for k in range(1,order+2):
            ans += ais[k-1,:]/k * (xis[-1]**k - xis[0]**k)
    return ans

使用了以下测试：

import numpy as np
from numpy import cos, sin , exp
import pyximport; pyximport.install()
import integratev
from subprocess import Popen
def func(x):
    return np.array([x**2, x**3, cos(x), sin(x), exp(x)])

if __name__ == '__main__':
    xs = np.linspace(0.,20.,33)
    print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.)+1, exp(20.)-1])
    ans =  integratev.trapzv(func,xs,5)
    print 'trapzv:', ans
    ans =  integratev.poly(func,xs,5,2)
    print 'poly:', ans

给予：

exact: [  2.66666667e+03   4.00000000e+04   9.12945251e-01   5.91917938e-01 4.85165194e+08]
trapzv: [  2.66796875e+03   4.00390625e+04   8.83031547e-01   5.72522998e-01 5.00856448e+08]
poly: [  2.66666675e+03   4.00000000e+04   9.13748980e-01   5.92435718e-01 4.85562144e+08]

多边形可以是任何顺序，这可能会在大多数情况下提供更好的结果......

于 2013-07-08T09:04:21.020 回答