Построение многомерной матрицы дифференцирования - PullRequest
Новая
       104

Построение многомерной матрицы дифференцирования

1 голос
/ 10 ноября 2019

Я пытался построить матрицу D ij , определенную как

definition image

Iхотите построить его для точек, расположенных в x i = -cos [ π (2 i + 1) / (2 N )] на интервале [-1,1], чтобы последовательно получить производные функции. Хотя у меня проблемы с построением дифференцирующей матрицы D ij .

Я написал скрипт на 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)

Таким образом, возвращается ошибка: 

IndexError: too many indices for array

Есть ли способ более эффективно построить это и успешновычислить это по всем к? Целью будет умножение этой матрицы на функцию и суммирование по j, чтобы получить производную первого порядка данной функции.

Ответы [ 3 ]

1 голос
/ 11 ноября 2019

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))
0 голосов
/ 12 ноября 2019

enter image description here

может быть реализовано как 

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]:
0 голосов
/ 11 ноября 2019

Посмотрите на свой x calc отдельно 

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 не выглядят интересными, потому что они всего лишь cos из нечетных целых кратных pi.

Обратите внимание, как я использую np.arange вместо циклана range.

...