Вопрос по numpy, scipy, python – Вычислить расхождение векторного поля с помощью Python

8

Есть ли функция, которая может быть использована для расчета расходимости векторного поля? (вMATLAB) Я ожидал бы, что это существует в numpy / scipy, но я не могу найти это, используя Google.

Мне нужно рассчитатьdiv[A * grad(F)], где

F = np.array([[1,2,3,4],[5,6,7,8]]) # (2D numpy ndarray)

A = np.array([[1,2,3,4],[1,2,3,4]]) # (2D numpy ndarray)

такgrad(F) это список 2Dndarrays

Я знаю, что могу рассчитать расхождение какэтот но не хочу изобретать велосипед. (Я бы также ожидал что-то более оптимизированное) У кого-нибудь есть предложения?

@mgilson Да, массивы расположены на одинаковом расстоянии. Мне нужна двойная точность. nyvltak
Какая точность заказа вам нужна? ваши массивы одинаково разнесены? mgilson
@ZagorulkinDmitry, Jensen & # x2013; Дивергенция Шеннона - это нечто совершенно другое nyvltak
@nyvltak - Не точность, порядок. Как вO(h) или жеO(h**2)и какой интервал? ... mgilson

Ваш Ответ

8   ответов
11

вышеприведенные функции не вычисляют расходимость векторного поля. они суммируют производные скалярного поля A:

result = dA/dx + dA/dy

в отличие от векторного поля (с трехмерным примером):

result = sum dAi/dxi = dAx/dx + dAy/dy + dAz/dz

Проголосуй за всех! Это математически просто неправильно.

Ура!

Как-то странно, как это внизу страницы. Другие ответы действительно математически неверны
2

http://www.scipy.org/Topical_Software#head-85e01502b533f2477ab8c643b38ee92706a377bb

Даже если он не имеет расходимости, упакованной для вас вручную, расхождение довольно простое, и производные инструменты, которые они дают вам в scipy (те, что указаны выше), дают вам около 90% кода, предварительно упакованного в хороший, эффективный путь.

1

ответ заключается в том, что в numpy нет встроенной функции дивергенции. Следовательноbest Метод вычисления расхождения заключается в суммировании компонентов вектора градиента, т.е.calculate расхождение.

1

но не Numpy. Это то, что, возможно, стоит внести в pylab, чтобы создать жизнеспособную альтернативу matlab с открытым исходным кодом.

http://wiki.scipy.org/PyLab

Изменить: теперь называетсяhttp://www.scipy.org/stackspec.html

9

def divergence(F):
    """ compute the divergence of n-D scalar field `F` """
    return reduce(np.add,np.gradient(F))

Timeit:

F = np.random.rand(100,100)
timeit reduce(np.add,np.gradient(F))
# 1000 loops, best of 3: 318 us per loop

timeit np.sum(np.gradient(F),axis=0)
# 100 loops, best of 3: 2.27 ms per loop

Примерно в 7 раз быстрее: sum неявно построить 3d-массив из списка полей градиента, которые возвращаютсяnp.gradient, Этого избегают, используяreduce

Теперь, в вашем вопросе, что вы подразумеваете подdiv[A * grad(F)]?

about A * grad(F): A is a 2d array, and grad(f) is a list of 2d arrays. So I considered it means to multiply each gradient field by A. about applying divergence to the (scaled by A) gradient field is unclear. By definition, div(F) = d(F)/dx + d(F)/dy + .... I guess this is just an error of formulation.

За1, умножая суммированные элементыBi тем же факторомA можно факторизовать:

Sum(A*Bi) = A*Sum(Bi)

Таким образом, вы можете получить этот взвешенный градиент просто с помощью:A*divergence(F)

Если & # x300;A вместо этого список факторов, один для каждого измерения, тогда решение будет:

def weighted_divergence(W,F):
    """
    Return the divergence of n-D array `F` with gradient weighted by `W`

    ̀`W` is a list of factors for each dimension of F: the gradient of `F` over
    the `i`th dimension is multiplied by `W[i]`. Each `W[i]` can be a scalar
    or an array with same (or broadcastable) shape as `F`.
    """
    wGrad = return map(np.multiply, W, np.gradient(F))
    return reduce(np.add,wGrad)

result = weighted_divergence(A,F)
Разве не расходимость векторного поля F = d (Fx) / dx + d (Fy) / dy + ...? Правильная формула будет что-то вродеnp.ufunc.reduce(np.add, [np.gradient(F[i], axis=i) for i in range(len(F))])
Это все зависит от типа данных в F. Вопрос неясен. У меня есть опыт обработки изображений, поэтому я считаю, что F - это изображение в формате nD, и поэтому градиент является производнойs отaxis х то у (и больше, если есть). Который я суммирую. Если я правильно понимаю, если F - это n * m 2D последовательность из n векторов, которые m-мерны, то я полагаю, что ваша формулировка верна. Однако я бы не понял, если F больше 2d
2

что Даниэль изменил, является правильным ответом, позвольте мне подробнее объяснить самоопределение функций:

Функция np.gradient () определяется как: np.gradient (f) = df / dx, df / dy, df / dz + ...

но нам нужно определить func расхождение как: расхождение (f) = dfx / dx + dfy / dy + dfz / dz + ... = np.gradient (fx) + np.gradient (fy) + np.gradient (fz) + ...

Давайте проверим, сравните спример расхождения в матлабе

import numpy as np
import matplotlib.pyplot as plt

NY = 50
ymin = -2.
ymax = 2.
dy = (ymax -ymin )/(NY-1.)

NX = NY
xmin = -2.
xmax = 2.
dx = (xmax -xmin)/(NX-1.)

def divergence(f):
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

y = np.array([ ymin + float(i)*dy for i in range(NY)])
x = np.array([ xmin + float(i)*dx for i in range(NX)])

x, y = np.meshgrid( x, y, indexing = 'ij', sparse = False)

Fx  = np.cos(x + 2*y)
Fy  = np.sin(x - 2*y)

F = [Fx, Fy]
g = divergence(F)

plt.pcolormesh(x, y, g)
plt.colorbar()
plt.savefig( 'Div' + str(NY) +'.png', format = 'png')
plt.show()

enter image description here

3

но изменено для правильной расходимости формулы векторного поля

def divergence(f):
    """
    Computes the divergence of the vector field f, corresponding to dFx/dx + dFy/dy + ...
    :param f: List of ndarrays, where every item of the list is one dimension of the vector field
    :return: Single ndarray of the same shape as each of the items in f, which corresponds to a scalar field
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

Документация Matlab использует эту точную формулу (прокрутите вниз до Расхождения векторного поля)

9
import numpy as np

def divergence(field):
    "return the divergence of a n-D field"
    return np.sum(np.gradient(field),axis=0)

Похожие вопросы