您的add适用于此计算，但必须重复应用，以便传播非零值。我不明白您的计算是如何生成[ 1. 1. 2. 1. 3. 1. 4. 1. 5. 1.]的。

In [619]: fib=np.zeros(10,int);fib[:2]=1
In [620]: for _ in range(10):
     ...:     np.add(fib[:-2], fib[1:-1], out=fib[2:])
     ...:     print(fib)
     ...: 
[1 1 2 1 0 0 0 0 0 0]   #**
[1 1 2 3 3 1 0 0 0 0]
[1 1 2 3 5 6 4 1 0 0]
[ 1  1  2  3  5  8 11 10  5  1]
[ 1  1  2  3  5  8 13 19 21 15]
 ...

(编辑-请注意，第一个np.add的作用就像它是完全缓冲的一样。将**处的结果与对象数组和浮点数组进行比较。我在Py3上使用的是1.13.1版。)
或者突出每个阶段的好的价值

In [623]: fib=np.zeros(20,int);fib[:2]=1
In [624]: for i in range(10):
     ...:     np.add(fib[:-2], fib[1:-1], out=fib[2:])
     ...:     print(fib[:(i+2)])
     ...: 
[1 1]
[1 1 2]
[1 1 2 3]
[1 1 2 3 5]
[1 1 2 3 5 8]
[ 1  1  2  3  5  8 13]
[ 1  1  2  3  5  8 13 21]
[ 1  1  2  3  5  8 13 21 34]
[ 1  1  2  3  5  8 13 21 34 55]
[ 1  1  2  3  5  8 13 21 34 55 89]

fib[2:] = fib[:-2]+fib[1:-1]做同样的事情。
正如ufunc.at的文档中所讨论的，像np.add这样的操作使用缓冲，即使使用out参数也是如此。因此，尽管它们在C级别代码中迭代，但它们不会累加；也就是说，ith步骤的结果不会在i+1步骤中使用。
add.at可用于执行无缓冲的a[i] += b。当索引包含DUPICATE时，这很方便。但它不允许从更改后的a值反馈到b。所以它在这里没有用处。
也是add.accumulate(也称为np.cumsum)，对于某些迭代定义很方便。但在一般情况下很难适用。
cumprod(乘法累加)可以使用qmatrix方法
Numpy's matrix_power function giving wrong results for large exponents
使用np.matrix定义qmatrix，因此*乘法表示矩阵乘积(而不是元素乘积)：

In [647]: qmatrix = numpy.matrix([[1, 1], [1, 0]])

用指向该矩阵的副本(实际上是指针)填充对象矩阵。

In [648]: M = np.empty(10, object)
In [649]: M[:] = [qmatrix for _ in range(10)]

应用cumprod，并从每个矩阵中提取一个元素。

In [650]: m1=np.cumprod(M)
In [651]: m1
Out[651]: 
array([matrix([[1, 1],
        [1, 0]]), matrix([[2, 1],
        [1, 1]]),
       matrix([[3, 2],
        [2, 1]]), matrix([[5, 3],
        [3, 2]]),
       matrix([[8, 5],
        [5, 3]]),
       matrix([[13,  8],
        [ 8,  5]]),
       matrix([[21, 13],
        [13,  8]]),
       matrix([[34, 21],
        [21, 13]]),
       matrix([[55, 34],
        [34, 21]]),
       matrix([[89, 55],
        [55, 34]])], dtype=object)
In [652]: [x[0,1] for x in m1]
Out[652]: [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]

链接答案使用numpy.linalg.matrix_power，它也做重复的矩阵乘积。对于单个电源，matrix_power更快，但对于整个电源范围，cumprod更快。

以下是使用Scipy过滤器的一种相当有效的方法：

import numpy as np
from scipy import signal

def scipy_fib(n):
    x = np.zeros(n, dtype=np.uint8)
    x[0] = 1

    sos = signal.tf2sos([1], [1, -1, -1])
    return signal.sosfilt(sos, x)

(请注意，我们不能使用lfilter，因为在信号处理方面，lfilter是一个不稳定的过滤器，会导致Python解释器崩溃，因此我们必须转换为二阶部分形式。)
过滤方法的问题是，我不确定它是否可以推广到解决一个ODE的实际问题。
另一种解决方案是简单地用Python编写循环，并使用JIT compiler加速它：

import numba as nb    

@nb.jit(nopython=True)
def numba_fib(n):
    y = np.empty(n)
    y[:2] = 1

    for i in range(2, n):
        y[i] = y[i-1] + y[i-2]

    return y

JIT方法的优点是您可以实现所有花哨的东西，但如果您坚持使用简单的循环和分支，而不调用任何(非JIT)的Python函数，那么它工作得最好。
计时结果：

# Baseline without JIT:

%timeit numba_fib(10000)  # 100 loops, best of 3: 5.47 ms per loop

%timeit scipy_fib(10000)  # 1000 loops, best of 3: 719 µs per loop
%timeit numba_fib(10000)  # 10000 loops, best of 3: 33.8 µs per loop

# For fun, this is how Paul Panzer's solve_banded approach compares:

%timeit banded_fib(10000)  # 1000 loops, best of 3: 1.33 ms per loop

Scipy Filter方法比纯Python循环快，但比JITted循环慢。我猜Filter函数是相对通用的，它做了我们在这里不需要的事情，而JITted循环被编译成一个非常小的循环。

一个可能不是超级有效但有趣的解决方案是像这样滥用scipy.linalg.solve_banded

import numpy as np
from scipy import linalg

N = 50
a = np.zeros((N,)) + np.array([[1, -1, -1]]).T
b = np.zeros((N,))
b[0] = 1
linalg.solve_banded((2, 0), a, b)

# array([  1.00000000e+00,   1.00000000e+00,   2.00000000e+00,

# 3.00000000e+00,   5.00000000e+00,   8.00000000e+00,

# 1.30000000e+01,   2.10000000e+01,   3.40000000e+01,

# 5.50000000e+01,   8.90000000e+01,   1.44000000e+02,

# 2.33000000e+02,   3.77000000e+02,   6.10000000e+02,

# 9.87000000e+02,   1.59700000e+03,   2.58400000e+03,

# 4.18100000e+03,   6.76500000e+03,   1.09460000e+04,

# 1.77110000e+04,   2.86570000e+04,   4.63680000e+04,

# 7.50250000e+04,   1.21393000e+05,   1.96418000e+05,

# 3.17811000e+05,   5.14229000e+05,   8.32040000e+05,

# 1.34626900e+06,   2.17830900e+06,   3.52457800e+06,

# 5.70288700e+06,   9.22746500e+06,   1.49303520e+07,

# 2.41578170e+07,   3.90881690e+07,   6.32459860e+07,

# 1.02334155e+08,   1.65580141e+08,   2.67914296e+08,

# 4.33494437e+08,   7.01408733e+08,   1.13490317e+09,

# 1.83631190e+09,   2.97121507e+09,   4.80752698e+09,

# 7.77874205e+09,   1.25862690e+10])

它的工作原理是我们写成F_0，F_1，F_2…作为一个向量F和恒等式-F_{i-1}-F_i+F{i+1}=0作为一个很好的带状矩阵，然后求F。请注意，这可以适用于其他类似的递推。

由于线性递归有解析解(这里是斐波那契)，另一种快速的方法是：

import math

import numpy as np

def fib_scipy2(N):
    r5=math.sqrt(5)
    phi=(1+r5)/2
    a= (-phi*np.ones(N)).cumprod()
    return (np.abs(a)-1/a)/r5

运行：

In [413]: fib_scipy2(8)
Out[413]: array([  1.,   1.,   2.,   3.,   5.,   8.,  13.,  21.])
In [414]: %timeit fib_scipy2(10**4)
103 µs ± 888 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

您可以将其调整为任何线性方程。