Численно решаем волновое уравнение разностной схемой15.03.2023 13:31

Для меня уравнения в частных производных — это очень красивая история из студенчества. Почему? Это невероятно красиво. Но что особенно стало для меня захватывающим, так это то, что дифуры в широком смысле прикладной математики — это тот самый пример, когда математика и компьютер используются вместе, чтобы представить некоторую компьютерную модель вполне реальных процессов. Как вы уже, наверное, догадались, речь пойдёт про то, как вообще можно попробовать решать дифференциальные уравнения в частных производных на компьютере. Мы попробуем это сделать на примере волнового уравнения и с использованием уже ставших привычными python, scipy и numpy. Если вы примерно помните математику, но панически боялись дифуров или они просто как-то обошли вас стороной, то добро пожаловать.

Волновое уравнение

Если просто, то дифференциальное уравнение — это некоторое уравнение на функцию, в которое входят также её производные. Много уже написано про обыкновенные дифференциальные уравнения, где неизвестная функция зависит от одного аргумента (например, см. здесь и здесь). Когда мы говорим об уравнениях в частных производных, то в них, соответственно, входят частные производные, также не обойдённые пользователями Хабра (зрелищно тут и тут). Вполне естественно, что решение таких уравнений не всегда существует, а если существует, то не обязательно единственное, а для единственности нам нужны некоторые дополнительные условия.

В данном случае мы рассматриваем волновое уравнение для функции u(t,x)

$\partial^2_t u = c^2 \partial^2_x u, \quad u: [0,T] \times [a,b] \to \mathbb{R}.$ Почему именно такое уравнение?

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

$m\ddot{u}_i = k (u_{i+1}-u_i) + k (u_{i-1}-u_{i}),$

где в правой части мы записали силу упругости в соответствии с законом Гука. Теперь заметим, что коэффициент жёсткости всей линии грузиков (это же система пружин) равен K=k/N , а масса грузика равна m=M/N , за обозначена масса всей системы. Подставим и немного перепишем:

$\ddot{u}_i = \frac{KN^2}{M} (u_{i+1}-2u_i + u_{i-1}).$

Что ещё постоянно? Длина между закреплёнными концами, так что горизонтальное расстояние между соседями равно h=L/N . Давайте теперь умножим и разделим правую часть на L^2 и…

$\ddot{u}_i = \frac{KL^2}{M} \frac{1}{h^2}(u_{i+1}-2u_i + u_{i-1}).$

Знакомое выражение справа? Да, теперь можем перейти к пределу при $h \to 0$ и получим уже известное

Дополнительно нам понадобятся начальные условия

$u(0,x)=u_0(x), ~ \partial_t u\vert_{(0,x)}=p(x)$

и в качестве краевых условий рассмотрим краевые условия Дирихле (функция как бы закреплена в нуле на краях)

$\forall t \in \mathbb{R}_+ ~ u(t,a)=0,~u(t,b)=0.$

Решением такой задачи будет динамика профиля натянутой cтруны, который изначально выглядел как u_0(x) , но под действием сил упругости изменяется и в момент выглядит как u(t,x) .

Конечно-разностная схема

Вообще, решение в данном случае даётся, например, формулой Д'Аламбера, но предположим, что мы не изучали уравнения в частных производных и сейчас этого не знаем. Какой разумный рецепт может прийти в голову? Вроде мы помним как примерно оценить производную со времён матана-1…

Почему бы и нет, но подойдём к вопросу строго. По сути мы пойдём в обратную сторону, пытаясь вообразить струну набором грузиков, соединённых пружиной. Давайте поделим струну на N_x равных сегментов и представим её как набор точек.Аналогичным образом поступим со временем, поделив интересующий нас отрезок [0,T] на N_t равных кусков; итого имеем N_x+1 точек x_j по пространству, N_t+1 точек t_i по времени, всего (N_x+1)(N_t+1) различных точек, образующих сетку $G_{N_x,N_t}$ . Наша задача будет состоять в том, чтобы построить такую дискретную (как ещё говорят, сеточную) функцию $\tilde{u}(t_i,x_j)$ , которая бы хорошо приближала настоящую u(t,x) . Конкретнее, нам бы хотелось, чтобы

$\max_{t_i,x_j \in G_{N_x,N_t}} \left\vert u(t_i,x_j) - \tilde{u}(t_i,x_j) \right\vert \to 0 \text{ при } N_x, N_t \to \infty.$

Расстояние между узлами дискретной сетки по обозначим h_x=(b-a)/N_x , а h_t=T/N_t .Так, при более густой дискретизации приближённое решение будет ближе к истинному.

Теперь разберёмся с дифференциальным уравнением

$\partial^2_t u = c^2 \partial^2_x u.$

Если пронумеруем точки по каждой координате от нуля, то у нас получатся индексы i=0,..,N_t и j=0,..,N_x . Заменим уравнение на похожее для нашей дискретной функции, оценив производную в точке (t_i,x_j) . Вспомним, что у нас есть начальное и краевые условия, поэтому новое уравнение мы записываем не во всех точках:

$\forall i\geq 1, 1\leq j \leq N_x-1 \\ \frac{\tilde{u}(t_{i+1},x_j) - 2\tilde{u}(t_i,x_j) + \tilde{u}(t_{i-1},x_j)}{h_t^2} = c^2 \frac{\tilde{u}(t_{i},x_{j+1}) - 2\tilde{u}(t_i,x_j) + \tilde{u}(t_{i},x_{j-1})}{h_x^2}.$

Мы знаем $\tilde{u}(t_0,x_j)$ т.к. в t_0=0 у нас есть начальное условие, таким же образом нам известны $\tilde{u}(t_i,x_{N_x})$ и $\tilde{u}(t_i,x_0)$ из краевых условий Дирихле. Есть, правда, небольшая проблема: нам неизвестно $\tilde{u}(t_1,x_j)$ , --, но её мы можем решить используя второе начальное условие. Помните, была первая производная в момент t=0 ? Вспомним старый добрый ряд Тэйлора и аппроксимируем u(t_1,x_j) :

$\tilde{u}(t_1,x_j) = \tilde{u}(0,x_j) + h_t \partial_t u \vert_{(0,x_j)}.$

Теперь остаётся только выразить $\tilde{u}(t_{i+1},x_{j})$ из уравнения выше и voilà, мы получаем явную схему (в том смысле, что есть простая итерационная формула $u(t_{i+1},\cdot)$ при известном $u(t_{i},\cdot)$ ). Наш алгоритм заключается в последовательном вычислении $\tilde{u}(t_2, \cdot),\tilde{u}(t_3, \cdot),...$ до самого последнего слоя по времени $\tilde{u}(t_{N_t}, \cdot)$ .

Подытожим, обозначив вектор $\tilde{u}_i:=\tilde{u}(t_i, \cdot)$ и введя матрицу размера $(N_x-1) \times (N_x-1)$ как матрицу с элементами

$\frac{2}{h_t^2} - \frac{2c^2}{h_x^2}$

на диагонали и c^2/h_x^2 над диагональю и под ней. Такая матрица называется трёхдиагональной. Наш алгоритм выглядит как-то так:

$\textbf{Вход:}~ u_0(x), \partial_t u \vert_{0,x}, N_x, N_t\\ \tilde{u}(0,x_j)\leftarrow u_0(x_j)\\ \tilde{u}(t_1,x_j) = \tilde{u}(0,x_j) + h_t \partial_t u \vert_{(0,x_j)} \;\\ for~ i \in \lbrace 2,...,N_t \rbrace ~~ \tilde{u}_i = h_t^2 C\tilde{u}_{i-1} - \tilde{u}_{i-2}$

Вычислений потребуется всего лишь порядка. Напишем небольшой класс, имплементирующий явную схему с таким решением. Для имплементации мы будем использовать scipy.sparse, который нам позволит задействовать разреженные матрицы, сэкономив много памяти и времени на умножения. Это будет особенно важно, если мы потом захотим решать волновое уравнение в размерности d>1» src=«https://habrastorage.org/getpro/habr/upload_files/57c/785/3dc/57c7853dc2642558a9712fa6f35af6ff.svg» />.</p>

<pre><code class= import numpy as np import scipy.sparse as sps import scipy.sparse.linalg as splin class ExplicitEulerSolver: def __init__(self, hs, cCoef): ''' Input float[] hs -- the steps of the discretization, [d+1], the first is the time ''' super(ExplicitEulerSolver,self).__init__() self.hs = hs self.cCoef = cCoef self.dim = hs.shape[0]-1 print(self.dim) def solveBVProblem1D(self, uInits, numSteps, rightPart=None, spaceGrid=None): ''' Solves Dirichlet Boundary-Value Problem Input float[] uInits -- two initial layers of shape [2,...dims...] int numSteps -- number of time steps funcHandler rightPart -- right part depending on time and space (t,x) float[][] spaceGrid -- space grid of points ''' #trim inits uSol = [ (uInits[0,...])[None,...], (uInits[1,...])[None,...] ] shape= np.array(list(uSol[0].shape)[1:])-2#-2 for dirichlet, trimmed shape for later restoration uSol[0] = np.take(uSol[0],np.arange(1,uSol[0].shape[1]-1),axis=1) uSol[1] = np.take(uSol[1],np.arange(1,uSol[1].shape[1]-1),axis=1) uSol[0] = np.reshape(uSol[0], newshape=(-1,), order='C')[None,:] uSol[1] = np.reshape(uSol[1], newshape=(-1,), order='C')[None,:] #assembling iteration matrix diags = self.hs[0]**2 *np.concatenate([self.cCoef**2 /(self.hs[1]**2) * np.ones([1,shape[0]]), (2/(self.hs[0]**2) - 2*self.cCoef**2 /(self.hs[1]**2) ) * np.ones([1,shape[0]]), self.cCoef**2 /(self.hs[1]**2) * np.ones([1,shape[0]])],axis=0) offsets = np.array([-1,0,1]) C = sps.spdiags(diags,offsets, m=len(diags[1]),n=len(diags[1]), format="csr") for i in np.arange(numSteps-1): if(rightPart is None): toAdd = C@uSol[-1][0,:] - uSol[-2][0,:] else: toAdd = C@uSol[-1][0,:] - uSol[-2][0,:] + np.reshape(rightPart(self.hs[0]*(i+2),spaceGrid),newshape=(-1,)) uSol.append(toAdd[None,:]) uSol=np.concatenate(uSol, axis=0) toShape= tuple([uSol.shape[0]] + [shape[j] for j in np.arange(len(shape))]) return np.pad(np.reshape(uSol, newshape=toShape), pad_width=[(0,0)]+[(1,1) for _ in np.arange(self.dim)])

Давайте проверим! Пусть a=0,b=2, T=4, c=2 , зададим начальные условия

$u_0(x)= \begin{cases} 2x, & x <0.5,\\ -\frac{2}{3}x + \frac{4}{3}, & x\geq 0.5, \end{cases} ~~\partial_t u\vert_{0,x}=0.$ Примерно соответствует тому, что мы просто оттянули струну

Примерно соответствует тому, что мы просто оттянули струну

Возьмём, к примеру, N_x=20,N_t=40 и нарисуем результат, используя следующий скрипт:

import numpy as np
import pickle

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation 

import LinearPDESolvers as linPDEs


import time
#solve wave eq

T=4
xLeft=0
xRight=2
c=2

#Nx = 510
#Nt = 1000
Nx=20
Nt=200
hx = (xRight-xLeft)/Nx
ht = T/Nt


solver = linPDEs.ExplicitEulerSolver(np.array([ht,hx]), cCoef=c)


xs = np.arange(xLeft,xRight+hx/2,hx)
ts = np.arange(0,T+ht/2,ht)

#assume u'(0,x)=0
def initFunc(xx):
    return ((2*xx)*(xx<0.5) + (xx>=0.5)*(-2/3*xx + 4/3))
def dtFunc(xx):
    return 0
uInits = np.zeros([2,len(xs)])
uInits[0,:] = initFunc(xs)
duts = dtFunc(xs)
uInits[1,:] = uInits[0,:] + duts*ht

#SOLVE!
res=solver.solveBVProblem1D(uInits, Nt)


# initializing a figure in 
# which the graph will be plotted
fig = plt.figure() 
   
# marking the x-axis and y-axis
axis = plt.axes(xlim =(xLeft, xRight), 
                ylim =(-1.5, 1.5)) 
  
# initializing a line variable
line, = axis.plot([], [], lw = 3) 

cross = axis.scatter([xLeft,xRight],[0,0], marker="x", s=90, color="red")

txt = axis.text(xLeft+1e-1,-1,"t="+str(ht*0))
# data which the line will 
# contain (x, y)
def init(): 
    line.set_data(xs, res[0,:])
    txt.set_text("t="+str(ht*0))
    return line,txt,
   
def animate(i):
    line.set_data(xs, res[i,:])
    txt.set_text("t={:.2f}".format(ht*i) )
    return line,txt,
   
anim = FuncAnimation(fig, animate, init_func = init,
                     frames = np.arange(0,res.shape[0],1), interval = ht*1000, blit = False, repeat=False)#ms scale

anim.save('explEuler.mp4',fps=30)

with open("./solveTempOut.pkl", "wb") as f:
    pickle.dump({"res": res},f)

Так стоп. Может, нужно попробовать более жесткую дискретизацию, сделаем, например, N_x=100 .

В реальности гитарные струны не разлетаются по Вселенной. В чём же дело? Может, слишком, грубая сетка, возьмём N_t=200 и…

вроде разумно выглядит. Но почему нам нужно так параноидально задавать шаг по времени, чтобы получить хоть что-то адекватное?

Устойчивость разностной схемы

Давайте вспомним наш итеративный метод; при известном $\tilde{u}_{i-1},\tilde{u}_{i-2}$

$\tilde{u}_i = h_t^2 C\tilde{u}_{i-1} - \tilde{u}_{i-2}.$

Наш взгляд неизбежно останавливается на том, что вообще-то $\tilde{u}_i$ линейно зависит от $\tilde{u}_{i-1},\tilde{u}_{i-2}$ . Действительно, в блочном виде

$\begin{bmatrix} \tilde{u}_{i}\\ \tilde{u}_{i-2} \end{bmatrix} = \begin{bmatrix} h_t^2 C & -I\\ 0 & I \end{bmatrix} \begin{bmatrix} \tilde{u}_{i-1}\\ \tilde{u}_{i-2} \end{bmatrix},$

где внизу мы просто записали очевидное равенство $\tilde{u}_{i-2}=\tilde{u}_{i-2}$ . Мы наблюдаем, что решение очень быстро растёт по модулю, значит, определённо, у этой матрицы посередине есть собственные вектора с собственным значением по модулю больше . Которые со временем просто взрываются. Это в чистом виде артефакт алгоритма, с физикой всё в порядке.

Но как проверить? Есть, например, теорема Гершгорина, однако она даёт обычно очень грубую оценку собственных чисел. Считать собственные числа? Но это дорого и не будем же мы для каждого N_x,N_t так проверять адекватность схемы. Хотелось бы получить простой критерий, в который мы могли бы подставить N_x,N_t и сказать сразу, будет ли алгоритм численно устойчивым.

Именно так и называется это свойство, которое упрощённо заключается в том, что при малых возмущениях $\tilde{u}_{i-1}+\delta_1,\tilde{u}_{i-2}+\delta_2$ полученный $\tilde{u}_i$ также мало отличается от $\tilde{u}_i$ , полученного без возмущений. Возмущения $\delta_1,\delta_2$ же раскладываются на собственные вектора? Если там будет даже с коэффициентом $10^{-12}$ собственный вектор с собственным значением, скажем, $\lambda=2$ , то уже где-то после шага i=25 гитарная струна улетит в небо.

Знаменитая теорема Лакса (оригинальная статья реально хорошая математически, посмотрите, если не видели) утверждает, что если численная схема устойчива, то есть, и аппроксимирует дифференциальное уравнение (в том смысле, что для полиномов от x,t степени не выше приближённая формула для уравнения точна), то численное решение будет равномерно по сетке сходиться к настоящему и ошибка будет порядка max(h_t,h_k)^k . Можно по-разному строить критерий устойчивости разностной схемы, но в любом случае нам надо как-то исследовать спектр матрицы

$\begin{bmatrix} h_t^2 C & -I\\ 0 & I \end{bmatrix}.$

Матрица h_t^2 C симметричная и трёхдиагональная с повторяющейся тройкой чисел. Поэтому, введя множество $\mathbb{X}= h_x\lbrace -N_x,..,N_x\rbrace$ , векторы $\tilde{u}$ рассмотреть как сигналы $\tilde{u}: \mathbb{X} \to \mathbb{R}$ (добавив нули из условия Дирихле по краям), а их умножение на матрицу h_t^2 C можно рассмотреть как свёртку c фильтром $f_C:\mathbb{X} \to \mathbb{R}$ , который определяется как

$f_C(y) = \begin{cases} \frac{c^2}{h_x^2}, ~~ y=\pm h_x,\\ \frac{2}{h_t^2}-\frac{2c^2}{h_x^2}, ~~ y=0,\\ 0, ~~ \text{в остальных случаях}. \end{cases}$

Так мы получим в другом виде

$\tilde{u}_{i} = h_t^2 f_C * \tilde{u}_{i-1} - \tilde{u}_{i-2}.$

Мы не добавили ничего нового, свёртка — всё равно линейное преобразование и у него есть собственные числа. Но зачем же нам ещё свёртки, если не для дискретного преобразования Фурье!

$\mathcal{F}\left[\tilde{u}_{i}\right](\xi) = h_t^2 \mathcal{F}\left[f_C\right](\xi) \mathcal{F}\left[ \tilde{u}_{i-1} \right](\xi) - \mathcal{F}\left[\tilde{u}_{i-2}\right](\xi).$

Поскольку мы можем взять обратное преобразование Фурье, а ещё прямое и обратное дискретные преобразования Фурье сохраняют расстояния, так что спектр матрицы по модулю меньше тогда и только тогда, когда решение полученного рекуррентного уравнения не будет расти по модулю для всех $\xi$ . Оно не такое сложное, как кажется, например, можем вычислить явно

$\mathcal{F}[f_C](\xi)= \frac{2}{h_t^2} - \frac{2c^2}{h_x^2} + \frac{c^2}{h_x^2} \left( e^{ih_x \xi} - e^{-ih_x \xi}\right) = \frac{2}{h_t^2} - \frac{2c^2}{h_x^2} + \frac{2c^2}{h_x^2} \cos(h_x \xi).$

Теперь вернёмся к записи итерации. Давайте для краткости положим

$F_{i}:=\mathcal{F}\left[\tilde{u}_{i}\right](\xi),$

тогда получим

$F_{i} = h_t^2 \mathcal{F}[f_C](\xi) F_{i-1} - F_{i-2}.$

Нам, таким образом, нужно, чтобы у характеристического уравнения

$\lambda^2 - h_t^2 \mathcal{F}[f_C](\xi) \lambda + 1=0$

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

По формулам Виета произведение корней равно единице;
Уравнение имеет вещественные коэффициенты, следовательно, у него может быть либо пара вещественных решений, либо пара комплексных, но сопряжённых.

Пункт 1 говорит, что если корни вещественные, то один из них точно больше единицы, а кратная единица для всех $\xi$ точно решением быть не может. С другой стороны, если имеем два комплексно сопряжённых корня, то они будут лежать ровно на единичной окружности. Наш случай! Это происходит тогда и только тогда, когда дискриминант отрицателен, то есть

$h_t^4 \mathcal{F}[f_C](\xi)^2 -4 <0.$

Дальше дело техники, которая нас приводит к

$-2 + \frac{c^2 h_t^2}{h_x^2} (1- \cos(h_x \xi)) < 0$

и далее для выполнения неравенства для всех $\xi$ можем потребовать

$\frac{ch_t}{h_x} <1.$

Число слева иногда называют числом Куранта; в результате мы получили достаточное условие устойчивости. Таким образом, если мы задаём в алгоритме такие шаги h_t, h_x , что

$\frac{ch_t}{h_x} < 1,$

полученная схема будет устойчивой и мы не увидим аномалий. Интересный вопрос состоит в том, что будет, если достигается равенство? В этом случае схема тоже будет устойчива, так как подставив

$\frac{ch_t}{h_x}=1$

в дискриминант и посчитав корни, мы увидим, что они по модулю меньше единицы.

Наконец, является ли это условие необходимым? Вдруг есть устойчивые схемы, для которых этот критерий не выполняется… Нет, такого не бывает. Если мы посмотрим, что происходит с дискриминантом при равенстве 1, то мы заметим, что у нас возникает точка $\xi$ , где дискриминант равен нулю и там кратный корень, по модулю меньше 1, который проблемой не является. Однако при числе Куранта больше 1 для некоторых $\xi$ возникнут некратные вещественные корни и это уже однозначно плохой случай. В терминах преобразования Фурье это будет означать, что какие-то частоты в ошибке будут экспоненциально усиливаться с каждой итерацией алгоритма.

Попробуйте схему с числом Куранта равным 1 и немного больше. Результат налицо.

А что насчёт уравнения в размерности d>1?

Физически не всё так очевидно, но мы получим аналогичную краевую задачу

$\partial_t^2 u = c^2 \Delta u ~~\forall x \in \Omega~ \forall t \in [0,T],\\ u(t,x)=0 ~~\forall x \in \partial \Omega ~\forall t \in (0,T]~,\\ u(0,x)=u_0(x),$

где $\Delta:= \sum_{i=1}^d \partial_{x^i}^2$ -- оператор Лапласа, а $\Omega, \partial \Omega$ -- интересующий нас регион пространства и его граница (например, куб [0,1]^d и его граница). Построенная нами схема Эйлера будет строиться точно так же, но интересно, что и анализ устойчивости мы можем провести аналогично, заменив одномерную свёртку на d-мерную, а одномерное преобразование Фурье на d-мерное. В полученной итеративной схеме на слое i-1 центральная точка $\tilde{u}(t_i,x^1_{j_1},..,x^d_{j_d})$ будет иметь вес

$2 - 2c^2h_t^2\sum_{i=1}^d\frac{1}{h_{x^j}^2},$

а точки, отходящие по пространству на 1 клетку, получат вес $c^2/h_{x^j}^2$ . И даже что-то типа числа Куранта мы тоже получим, выведя критерий устойчивости, так же записав условие для дискриминанта (не буду утомлять формулами, да и мне кажется, можно получше результат получить):

$\forall i \in \lbrace 1,..,d \rbrace ~~ \frac{dc^2h_t^2}{h_{x^i}^2} < 1,$

что в случае d=1 даёт то, что мы получили ранее.

А что дальше?

Если пойти немного дальше и подумать про конструкцию разностной схемы, то можно построить такие, которые будут устойчивы при любом выборе шага. Вместе с теоремой Лакса это дало бы нам возможность не быть такими параноидальными в плане выборе шага по времени. С другой стороны, построенная явная схема Эйлера — это всего лишь схема первого порядка и для того, чтобы увидеть (более красивую картинку, см. ниже) много деталей, шаг по времени нужно задавать достаточно большим, что будет приводить к большим вычислительным затратам. Возможно, в следующий раз построим что-нибудь покруче. Да и уравнение в размерности d>1 тоже интересно решить.

Думаете, слишком регулярно? Посмотрите на реальное видео: