构建多维微分矩阵
Constructing a Multidimensional Differentiation Matrix
我一直在尝试构造矩阵Dij，定义为
我想为位于 xi = -cos[ π[=35= 的点绘制它] (2 i + 1) / (2 N)] 在区间 [-1,1] 上对函数求导。我在构建微分矩阵 Dij 时遇到了问题。
我写了一个 python 脚本作为：
import numpy as np
N = 100
x = np.linspace(-1,1,N-1)
for i in range(0, N - 1):
x[i] = -np.cos(np.pi*(2*i + 1)/2*N)
def Dmatrix(x,N):
m_ij = np.zeros(3)
for k in range(len(x)):
for j in range(len(x)):
for i in range(len(x)):
m_ij[i,j,k] = -2/N*((k*np.sin(k*np.pi*(2*i + 1)/2*N(np.cos(k*np.pi*(2*j +1))/2*N)/(np.sin(np.pi*(2*i + 1)/2*N)))
return m_ij
xx = Dmatrix(x,N)
因此 returns 错误：
IndexError: too many indices for array
有没有一种方法可以更有效地构造它并成功计算所有 k ？
目标是将这个矩阵乘以一个函数并对 j 求和以获得给定函数的一阶导数。
m_ij = np.zeros(3)不是三维数组，而是长度为3的一维数组。
In [1]: import numpy as np
In [2]: m_ij = np.zeros(3)
In [3]: print(m_ij)
[0. 0. 0.]
我怀疑你想要（作为一个简单的修复）
len_x = len(x)
m_ij = np.zeros((len_x, len_x, len_x))
单独查看您的 x 计算结果
In [418]: N = 10
...: x = np.linspace(-1,1,N-1)
...: y = np.zeros(N)
...: for i in range(N):
...: y[i] = -np.cos(np.pi*(2*i + 1)/2*N)
...:
In [419]: x
Out[419]: array([-1. , -0.75, -0.5 , -0.25, 0. , 0.25, 0.5 , 0.75, 1. ])
In [420]: y
Out[420]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
In [421]: (2*np.arange(N)+1)
Out[421]: array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19])
In [422]: (2*np.arange(N)+1)/2*N
Out[422]: array([ 5., 15., 25., 35., 45., 55., 65., 75., 85., 95.])
我把 x 和 y 分开了，否则创建 x 然后覆盖它就没有任何意义了。
y 值看起来并不有趣，因为它们都是 pi.
的奇数倍数 cos
请注意我是如何使用 np.arange 而不是在 range 上循环的。
可以实现为
def D(N):
from numpy import zeros, pi, sin, cos
D = zeros((N, N))
for i in range(N):
for j in range(N):
for k in range(N):
D[i,j] -= k*sin(k*pi*(i+i+1)/2/N)*cos(k*pi*(j+j+1)/2/N)
D[i,j] /= sin(pi*(i+i+1)/2/N)
return D*2/N
向量化内循环会很方便。
在第二次尝试中，所有过程都可以使用 np.einsum 进行矢量化（最后我也有一些时间，einsum 版本当然比三重循环快得多）：
In [1]: from numpy import set_printoptions ; set_printoptions(linewidth=120)
In [2]: def D(N):
...: from numpy import zeros, pi, sin, cos
...: D = zeros((N, N))
...: for i in range(N):
...: for j in range(N):
...: for k in range(N):
...: D[i,j] -= k * sin(k*pi*(2*i+1)/2/N) * cos(k*pi*(2*j+1)/2/N)
...: D[i,j] /= sin(pi*(2*i+1)/2/N)
...: return D*2/N
In [3]: def E(N):
...: from numpy import arange, cos, einsum, outer, pi, sin
...: i = j = k = arange(N)
...: s_i = sin((2*i+1)*pi/2/N)
...: s_ki = sin(outer(k,(2*i+1)*pi/2/N))
...: c_kj = cos(outer(k,(2*j+1)*pi/2/N))
...: return -2/N*einsum('k, ki, kj -> ij', k, s_ki, c_kj) / s_i[:,None]
In [4]: for N in (3,4,5):
...: print(D(N)) ; print(E(N)) ; print('==========')
...:
[[-1.73205081e+00 2.30940108e+00 -5.77350269e-01]
[-5.77350269e-01 1.22464680e-16 5.77350269e-01]
[ 5.77350269e-01 -2.30940108e+00 1.73205081e+00]]
[[-1.73205081e+00 2.30940108e+00 -5.77350269e-01]
[-5.77350269e-01 1.22464680e-16 5.77350269e-01]
[ 5.77350269e-01 -2.30940108e+00 1.73205081e+00]]
==========
[[-3.15432203 4.46088499 -1.84775907 0.5411961 ]
[-0.76536686 -0.22417076 1.30656296 -0.31702534]
[ 0.31702534 -1.30656296 0.22417076 0.76536686]
[-0.5411961 1.84775907 -4.46088499 3.15432203]]
[[-3.15432203 4.46088499 -1.84775907 0.5411961 ]
[-0.76536686 -0.22417076 1.30656296 -0.31702534]
[ 0.31702534 -1.30656296 0.22417076 0.76536686]
[-0.5411961 1.84775907 -4.46088499 3.15432203]]
==========
[[-4.97979657e+00 7.20682930e+00 -3.40260323e+00 1.70130162e+00 -5.25731112e-01]
[-1.05146222e+00 -4.49027977e-01 2.10292445e+00 -8.50650808e-01 2.48216561e-01]
[ 3.24919696e-01 -1.37638192e+00 2.44929360e-16 1.37638192e+00 -3.24919696e-01]
[-2.48216561e-01 8.50650808e-01 -2.10292445e+00 4.49027977e-01 1.05146222e+00]
[ 5.25731112e-01 -1.70130162e+00 3.40260323e+00 -7.20682930e+00 4.97979657e+00]]
[[-4.97979657e+00 7.20682930e+00 -3.40260323e+00 1.70130162e+00 -5.25731112e-01]
[-1.05146222e+00 -4.49027977e-01 2.10292445e+00 -8.50650808e-01 2.48216561e-01]
[ 3.24919696e-01 -1.37638192e+00 2.44929360e-16 1.37638192e+00 -3.24919696e-01]
[-2.48216561e-01 8.50650808e-01 -2.10292445e+00 4.49027977e-01 1.05146222e+00]
[ 5.25731112e-01 -1.70130162e+00 3.40260323e+00 -7.20682930e+00 4.97979657e+00]]
==========
In [5]: %timeit D(20)
36 ms ± 277 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [6]: %timeit E(20)
146 µs ± 777 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [7]: %timeit D(100)
4.35 s ± 30.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [8]: %timeit E(100)
7.7 ms ± 2.82 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [9]:
我一直在尝试构造矩阵Dij，定义为
我想为位于 xi = -cos[ π[=35= 的点绘制它] (2 i + 1) / (2 N)] 在区间 [-1,1] 上对函数求导。我在构建微分矩阵 Dij 时遇到了问题。
我写了一个 python 脚本作为：
import numpy as np
N = 100
x = np.linspace(-1,1,N-1)
for i in range(0, N - 1):
x[i] = -np.cos(np.pi*(2*i + 1)/2*N)
def Dmatrix(x,N):
m_ij = np.zeros(3)
for k in range(len(x)):
for j in range(len(x)):
for i in range(len(x)):
m_ij[i,j,k] = -2/N*((k*np.sin(k*np.pi*(2*i + 1)/2*N(np.cos(k*np.pi*(2*j +1))/2*N)/(np.sin(np.pi*(2*i + 1)/2*N)))
return m_ij
xx = Dmatrix(x,N)
因此 returns 错误：
IndexError: too many indices for array
有没有一种方法可以更有效地构造它并成功计算所有 k ？ 目标是将这个矩阵乘以一个函数并对 j 求和以获得给定函数的一阶导数。
m_ij = np.zeros(3)不是三维数组，而是长度为3的一维数组。
In [1]: import numpy as np
In [2]: m_ij = np.zeros(3)
In [3]: print(m_ij)
[0. 0. 0.]
我怀疑你想要（作为一个简单的修复）
len_x = len(x)
m_ij = np.zeros((len_x, len_x, len_x))
单独查看您的 x 计算结果
In [418]: N = 10
...: x = np.linspace(-1,1,N-1)
...: y = np.zeros(N)
...: for i in range(N):
...: y[i] = -np.cos(np.pi*(2*i + 1)/2*N)
...:
In [419]: x
Out[419]: array([-1. , -0.75, -0.5 , -0.25, 0. , 0.25, 0.5 , 0.75, 1. ])
In [420]: y
Out[420]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
In [421]: (2*np.arange(N)+1)
Out[421]: array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19])
In [422]: (2*np.arange(N)+1)/2*N
Out[422]: array([ 5., 15., 25., 35., 45., 55., 65., 75., 85., 95.])
我把 x 和 y 分开了，否则创建 x 然后覆盖它就没有任何意义了。
y 值看起来并不有趣，因为它们都是 pi.
cos
请注意我是如何使用 np.arange 而不是在 range 上循环的。
可以实现为
def D(N):
from numpy import zeros, pi, sin, cos
D = zeros((N, N))
for i in range(N):
for j in range(N):
for k in range(N):
D[i,j] -= k*sin(k*pi*(i+i+1)/2/N)*cos(k*pi*(j+j+1)/2/N)
D[i,j] /= sin(pi*(i+i+1)/2/N)
return D*2/N
向量化内循环会很方便。
在第二次尝试中，所有过程都可以使用 np.einsum 进行矢量化（最后我也有一些时间，einsum 版本当然比三重循环快得多）：
In [1]: from numpy import set_printoptions ; set_printoptions(linewidth=120)
In [2]: def D(N):
...: from numpy import zeros, pi, sin, cos
...: D = zeros((N, N))
...: for i in range(N):
...: for j in range(N):
...: for k in range(N):
...: D[i,j] -= k * sin(k*pi*(2*i+1)/2/N) * cos(k*pi*(2*j+1)/2/N)
...: D[i,j] /= sin(pi*(2*i+1)/2/N)
...: return D*2/N
In [3]: def E(N):
...: from numpy import arange, cos, einsum, outer, pi, sin
...: i = j = k = arange(N)
...: s_i = sin((2*i+1)*pi/2/N)
...: s_ki = sin(outer(k,(2*i+1)*pi/2/N))
...: c_kj = cos(outer(k,(2*j+1)*pi/2/N))
...: return -2/N*einsum('k, ki, kj -> ij', k, s_ki, c_kj) / s_i[:,None]
In [4]: for N in (3,4,5):
...: print(D(N)) ; print(E(N)) ; print('==========')
...:
[[-1.73205081e+00 2.30940108e+00 -5.77350269e-01]
[-5.77350269e-01 1.22464680e-16 5.77350269e-01]
[ 5.77350269e-01 -2.30940108e+00 1.73205081e+00]]
[[-1.73205081e+00 2.30940108e+00 -5.77350269e-01]
[-5.77350269e-01 1.22464680e-16 5.77350269e-01]
[ 5.77350269e-01 -2.30940108e+00 1.73205081e+00]]
==========
[[-3.15432203 4.46088499 -1.84775907 0.5411961 ]
[-0.76536686 -0.22417076 1.30656296 -0.31702534]
[ 0.31702534 -1.30656296 0.22417076 0.76536686]
[-0.5411961 1.84775907 -4.46088499 3.15432203]]
[[-3.15432203 4.46088499 -1.84775907 0.5411961 ]
[-0.76536686 -0.22417076 1.30656296 -0.31702534]
[ 0.31702534 -1.30656296 0.22417076 0.76536686]
[-0.5411961 1.84775907 -4.46088499 3.15432203]]
==========
[[-4.97979657e+00 7.20682930e+00 -3.40260323e+00 1.70130162e+00 -5.25731112e-01]
[-1.05146222e+00 -4.49027977e-01 2.10292445e+00 -8.50650808e-01 2.48216561e-01]
[ 3.24919696e-01 -1.37638192e+00 2.44929360e-16 1.37638192e+00 -3.24919696e-01]
[-2.48216561e-01 8.50650808e-01 -2.10292445e+00 4.49027977e-01 1.05146222e+00]
[ 5.25731112e-01 -1.70130162e+00 3.40260323e+00 -7.20682930e+00 4.97979657e+00]]
[[-4.97979657e+00 7.20682930e+00 -3.40260323e+00 1.70130162e+00 -5.25731112e-01]
[-1.05146222e+00 -4.49027977e-01 2.10292445e+00 -8.50650808e-01 2.48216561e-01]
[ 3.24919696e-01 -1.37638192e+00 2.44929360e-16 1.37638192e+00 -3.24919696e-01]
[-2.48216561e-01 8.50650808e-01 -2.10292445e+00 4.49027977e-01 1.05146222e+00]
[ 5.25731112e-01 -1.70130162e+00 3.40260323e+00 -7.20682930e+00 4.97979657e+00]]
==========
In [5]: %timeit D(20)
36 ms ± 277 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [6]: %timeit E(20)
146 µs ± 777 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [7]: %timeit D(100)
4.35 s ± 30.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [8]: %timeit E(100)
7.7 ms ± 2.82 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [9]: