Функция на сетке

Важная особенность np.array — возможность обходиться при вычислении значения функции без циклов.

import numpy as np

def f(x,y):
    return x**2+2*x+6*y

s = np.linspace(1,7,3)
p = np.linspace(5,9,3)

v = f (s,p) # Функция f была применена поэлементно 
# v == array([  33.,   66.,  117.])

S, P = np.meshgrid(s,p,indexing='ij')
v = f (S,P) # Функция f была применена для каждой пары из s и p
# array([[ 33.,  45.,  57.],
#       [ 54.,  66.,  78.],
#       [ 93., 105., 117.]])

В случае простой подстановки f (s,p) будет массивы будут подставлены поэлементно. В случае подстановки с использованием meshgrid из массивов будет построена сетка, в узлах которой будут пары элементов и функция будет рассчитана для каждого элемента сетки.

Дифференциальные операторы на сетке

Для расчетов производных на сетке будем использовать модуль findiff^pip (документация).

import numpy as np
from findiff import FinDiff

x = np.linspace(0, 10, 100)
dx = x[1] - x[0]
f = np.sin(x)

d_dx = FinDiff(0, dx, 1)
df_dx = d_dx(f)

Объект FinDiff реализует дифференциальный оператор, аргументы: 0 — по первой оси, dx — размер шага, 1 — первая производная.

df_dx — массив np.array со значениями производной f в узлах сетки. Модуль findiff не изменяет размер массива т.к. на краях применяется схема дифференциирование вперед или назад соответственно, а по остальному объему массива используется центральная схема.

Дифференциальные операторы модуля findiff

В 2D задаче

import numpy as np
from findiff import Gradient, Divergence, Laplacian
import matplotlib.pyplot as plt

#%%

borders = [-2*np.pi, 2*np.pi, -2*np.pi, 2*np.pi]

x = np.linspace(-2*np.pi, 2*np.pi, 30)
dx = x[1]-x[0]
y = np.linspace(-2*np.pi, 2*np.pi, 30)
dy = y[1]-y[0]

X, Y = np.meshgrid(x, y, indexing='ij')
f = Y*np.sin(X) + np.cos(Y)

#%%

plt.imshow(f.T,  origin='lower', extent=borders)
contours = plt.contour(X, Y, f,  origin='lower', colors = 'r')
plt.clabel(contours, colors = 'r')
plt.show()

#%%
grad = Gradient(h=[dx, dy])
grad_f = grad(f)

plt.imshow(f.T ,  origin='lower', extent=borders)
plt.colorbar()
plt.quiver(X,Y,grad_f[0],grad_f[1])
plt.show()

#%%

u = np.array((X**2-Y**2, Y**2-X))

div = Divergence(h=[dx, dy])
div_u = div(u)

plt.imshow(div_u.T,origin='lower', extent=borders)
plt.colorbar()
plt.quiver(X,Y,u[0],u[1])
plt.show()

Скалярная функция f и ее изолинии

Скалярная функция f и ее градиент

Векторная функция u и ее дивергенция

В 3D задаче

import numpy as np
from findiff import Gradient, Divergence, Laplacian, Curl

x = np.linspace(0, 10, 100)
dx = x[1]-x[0]
y = np.linspace(0, 10, 100)
dy = y[1]-y[0]
z = np.linspace(0, 10, 100)
dz = z[1]-z[0]

X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
f = np.sin(X) * np.cos(Y) * np.sin(Z)

grad = Gradient(h=[dx, dy, dz])
grad_f = grad(f)

laplace = Laplacian(h=[dx, dy, dz])
laplace_f = laplace(f)

g = np.array([f, 2*f, 3*f])

div = Divergence(h=[dx, dy, dz])
div_g = div(g)

curl = Curl(h=[dx, dy, dz])
curl_g = curl(g)