Эллиптический тип

• Обычно нет временной зависимости - несколько пространственных координат. Нет эволюции во времени - процесс стационарный.

• Решаются как краевые задачи при заданных краевых условиях.

• Решения гладкие, если позволяют краевые условия.

Примеры:

• Уравнение Пуассона (Лапласа при \(f=0\)):

\[ \sum_{i=1}^n \frac{\partial^2}{\partial x_i^2} u(x)=f(x) \]

Пример решения уравнения Пуассона в 2d. Прямой метод.

Мы хотим решить уравнение Пуассона:

\[\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x, y)\]

в \([0,1] \times [0,1]\) - квадратная область с граничными условиями

\[u(x,0) = x, \quad u(x,1) = x - 1, \quad u(0,y) = -y, \quad u(1,y) = 1 - y.\]

Примечание. Даже если область другая, можем произвести замену координат к квадрату. Или выпендриваться, о дискретизации нетривиальных граничных условий позже.

Дискретизуем производные. Возьмём формулы второго порядка аппроксимации по каждой из осей.

\[ \frac{u_{i-1, j}-2 u_{i, j}+u_{i+1, j}}{h^2}+\frac{u_{i, j-1}-2 u_{i, j}+u_{i, j+1}}{h^2}=f_{i j} \]

Отсюда получаем уравнение на центральную точку (т.к. знаем края логично двигаться в центр)

\[ u_{i, j}=\frac{1}{4}\left(u_{i+1, j}+u_{i-1, j}+u_{i, j+1}+u_{i, j-1}-h^2 f\left(x_i, y_j\right) \right) \]

Решаем наш «диффур» на сетке \((i, j)\). При этом знаем все граничные точки на сетке, и есть уравнение связи точек на «кресте» (т.н. шаблон крест).

Примечание. Зная 4 точки из «креста», можно найти пятую, выразив из уравнения. При другой дискретизации производных получили бы другие шаблоны.

Такая схема обладает 2-м порядком аппроксимации (достаточно просто разложить по Тейлору вокруг точек \(x_i\) и \(y_j\) и убедиться, что дифференциальное уравнение выполняется с точностью \(O(h^2)\)).

Для решения задачи (определения всех внутренних точек, зная только краевые) можно составить большую СЛАУ на все внутренние точки. Грубо говоря, будем идти нашим крестом по всем внутренним точкам начиная с нижнего ряда и до верхнего по порядку. В каждом случае, будем запоминать, какая функциональная связь между коэффициентами. Когда дойдём до верхнего ряда, пойдём назад и «расшифруем» все пройдённые точки.

Звучит сложно, да? Это называется прямой метод. Объясним его подробнее.

Хотим свести дифференциальную задачу к дискретной, заданной в тензорном (дискретизованном) виде (с учётом краевых условий):

\[ L_{i,j,k,l}\cdot u_{k,l}=b_{i,j} \]

Во всех «внутренних» точках квадрата оно имеет тот самый дискретизованный вид, который мы написали ранее:

\[ b_{i,j}=f_{i,j}, \quad i,j \in [1, \cdots,n-1] \]

\[ L_{i,j,k,l}\cdot u_{k,l}=\frac{u_{i-1, j}-2 u_{i, j}+u_{i+1, j}}{h^2}+\frac{u_{i, j-1}-2 u_{i, j}+u_{i, j+1}}{h^2}, \quad i,j \in [1, \cdots,n-1] \]

\[ L_{i,j,k,l}\cdot u_{k,l}=\frac{1}{h^2}\left(u_{i-1, j}+u_{i+1, j} +u_{i, j-1}+u_{i, j+1}\right)-\frac{4}{h^2}u_{i, j}, \quad i,j \in [1, \cdots,n-1] \]

А тогда просто вспомнив, как работает сумма, для внутренних \(i,j \in [1, \cdots,n-1]\) выполнено

\[\begin{split} \left\{\begin{array}{l} L_{i, j, i+1, j}=L_{i, j, i-1, j}=L_{i, j, i, j+1}=L_{i, j, i, j-1}=\frac{1}{h^2} \\ L_{i, j, i, j}=-\frac{4}{h^2} \\ L_{i, j, \text { other }}=0 \end{array}\right. \end{split}\]

Также в это тензорное уравнение мы хотим занести краевые условия. Достаточно следующего:

\[ b_{0,j}=\varphi_{\text{left}}(y_j) \]

\[ L_{0, j, k, l}\cdot u_{k, l} = u_{0, j} \]

что эквивалентно

\[\begin{split} \left\{\begin{array}{l} L_{0, j, 0, j}=1 \\ L_{0, j, \text{other}}=0 \end{array}\right. \end{split}\]

Для остальных трёх краевых условий аналогично.

Теперь, когда мы определили все элементы \(L_{i,j,k,l}\) и \(b_{i,j}\), остаётся только решить тензорное уравнение относительно \(u_{i,j}\):

\[ L_{i,j,k,l}\cdot u_{k,l}=b_{i,j} \]

В этом могут помочь методы numpy ;).

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt


def boundary(grid):
    x = np.linspace(0,1,len(grid))
    
    grid[0,:]  = -x
    grid[:,-1] = x-1
    grid[-1,:] = 1-x
    grid[:,0]  = x
    

def poisson_direct(gridsize, set_boundary):
    L = np.zeros(shape=(gridsize, gridsize, gridsize, gridsize),
                 dtype='d')
    
    b = np.zeros(shape=(gridsize, gridsize),
                 dtype='d')
    
    dx = 1.0 / (gridsize - 1)
    
    # discretized differential operator
    for i in range(1,gridsize-1):
        for j in range(1,gridsize-1):
            L[i,j,i-1,j] = L[i,j,i+1,j] = L[i,j,i,j-1] = L[i,j,i,j+1] = 1/dx**2
            L[i,j,i,j] = -4/dx**2
    
    # boundary conditions
    for i in range(0,gridsize):
        L[0,i,0,i] = L[-1,i,-1,i] = L[i,0,i,0] = L[i,-1,i,-1] = 1
    
    # set the boundary values on the right side
    set_boundary(b)
    
    return np.linalg.tensorsolve(L,b) 

def showsol(sol):
    plt.imshow(sol.T,cmap=mpl.cm.Blues,interpolation='none',origin='lower')

sol = poisson_direct(25, boundary) # Сетка из 25 точек
showsol(sol)

showsol(poisson_direct(50,boundary)) # Уменьшим сетку

Итерационный метод - сведение к уравнению теплопроводности

Вообще, прямой подход достаточно сложный - нужно обращать сильно разряженную матрицу. Можно решить эту же задачу итерационно.

На основе нашего уравнения Пуассона запишем параболическое уравнение диффузии (технически, оно уже от трёх переменных - см. увеличенную классификацию):

\[ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}-f(x,y) \]

Идея заключается в том, что на бесконечности по \(t\) эволюция системы должна остановится (система достигнет равновесия), а тогда левая часть занулится и полученное состояние будет решением уравнения Пуассона в чистом виде. Стартовать при этом «казалось бы» можно из любой начальной конфигурации, отвечающей граничным условиям. Естественно, скорость сходимости будет зависеть от близости к истинному решению.

Такое уже решать умеем - это же параболическое уравнение. Решим его явной схемой для простоты.

\[ u^{n+1}_{i,j} = \frac{1}{4}(u^{n}_{i+1,j} + u^{n}_{i-1,j} + u^{n}_{i,j+1} + u^{n}_{i,j-1})-\frac{h^2}{4}\cdot f_{i,j} \]

from numba import njit

@njit()
def boundary(grid):
    x = np.linspace(0,1,len(grid))
    
    grid[0,:]  = -x
    grid[:,-1] = x-1
    grid[-1,:] = 1-x
    grid[:,0]  = x


@njit()
def step(grid):
    newgrid = grid.copy()
    
    for i in range(1,newgrid.shape[0]-1):
        for j in range(1,newgrid.shape[1]-1):
            newgrid[i,j] = 0.25 * (newgrid[i,j+1] + newgrid[i,j-1] +
                                   newgrid[i+1,j] + newgrid[i-1,j])
    
    return newgrid


@njit() # Начальная случайная сетка
def initgrid(gridsize): 
    x = np.random.randn(gridsize,gridsize)
    boundary(x)
    
    return x

x = initgrid(100)

plt.figure(figsize=(12,12))

for i in range(1500):    
    if i % 100 == 0:
        plt.subplot(4,4,i//100+1)
        showsol(x)
        plt.title('iter = %s' % i)
    x = step(x)

plt.tight_layout()

# Посмотрим также на ошибку от временной итерации

x = initgrid(25)
sol = poisson_direct(25, boundary)

err = []

for i in range(1000):    
    err.append((i,np.linalg.norm(x - sol)))
            
    x = step(x)

err = np.array(err)
    
plt.semilogy(err[:,0],err[:,1])
plt.title('Error after iteration');