0

如何使用 python Numpy 矢量化解决第 2 页上的练习 4.5？

下载链接：

https://dl.dropbox.com/u/92795325/Python%20Scripting%20for%20Computational%20Scien%20-%20H.P.%20%20Langtangen.pdf

我用 Python 循环试过这个，但我需要矢量化版本。

from numpy import *
import time

def fun1(x):
       return 2*x+1

def integ(a,b,n):
       t0 = time.time()
       h = (b-a)/n
       a1 = (h/2)*fun1(a)
       b1 = (h/2)*fun1(b)
       c1 = 0
       for i in range(1,n,1):
              c1 = fun1((a+i*h))+c1
       t1 = time.time()
       return a1+b1+h*c1, t1-t0
4

1 回答 1

0

要“矢量化”使用numpy，这意味着不是像这样进行显式循环，

for i in range(1, n):
    c = c + f(i)

然后你应该把i它做成一个 numpy 数组，然后简单地取它的总和：

i = np.arange(1,n)
c = i.sum()

numpy 会自动为您进行矢量化。这更快的原因是由于各种原因，numpy 循环以比普通 python 循环更好的优化方式完成。一般来说，循环/数组越长，优势就越大。这是您实现的梯形积分：

import numpy as np

def f1(x):
    return 2*x + 1

# Here's your original function modified just a little bit:
def integ(f,a,b,n):
    h = (b-a)/n
    a1 = (h/2)*f(a)
    b1 = (h/2)*f(b)
    c1 = 0
    for i in range(1,n,1):
        c1 = f((a+i*h))+c1
    return a1 + b1 + h*c1

# Here's the 'vectorized' function:
def vinteg(f, a, b, n):
    h = (b-a) / n
    ab = 0.5 * h * (f(a)+f(b)) #only divide h/2 once

    # use numpy to make `i` a 1d array:
    i = np.arange(1, n) 
    # then, passing a numpy array to `f()` means that `f` returns an array
    c = f(a + h*i) # now c is a numpy array

    return ab + c.sum() # ab + np.sum(c) is equivalent

在这里，我将把我命名的内容导入tmp.py到ipython会话中，以便比使用更容易计时time.time：

import trap
f = trap.f1
a = 0
b = 100
n = 1000

timeit trap.integ(f, a, b, n)
#1000 loops, best of 3: 378 us per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 51.6 us per loop

哇，快七倍。

看看它是否对较小的人有很大帮助n

n = 10

timeit trap.integ(f, a, b, n)
#100000 loops, best of 3: 6 us per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 43.4 us per loop

不，小循环要慢得多！非常大n呢？

n = 10000

timeit trap.integ(f, a, b, n)
#100 loops, best of 3: 3.69 ms per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 111 us per loop

快三十倍！

于 2013-04-06T22:20:09.890 回答