Применение вариационных принципов для построения разностных схем

Имеем интегро - дифференциальное уравнение для определения $$\theta (s, t)$$, причем как строить его разностную аппроксимацию — непонятно.

Отметим, что $$G(s, \sigma )$$ есть функция Грина для задачи

w'' = - g(s),     w'(0) = w(1) = 0.

(рис 9.2)

Теперь введем дискретный аналог лагранжиана Lh. Для этого разобьем стержень на n отрезков одинаковой длины $$\Delta l$$ (и одинаковой массы). Каждый отрезок характеризуется углом наклона $$\theta _{k}$$. Тогда (см. рис. 9.2)

$$\begin{gather*} x_j = \sum\limits_{k = 1}^{j}{\Delta l \cos \theta_k } = \Delta l \sum\limits_{k = 1}^{j}{\cos \theta_k } ; \\ y_j = \Delta l \sum\limits_{k = 1}^{j}{{\sin}\theta_k } ; \\ T_{h} = \frac{\Delta l}{2} \sum\limits_{j = 1}^{n}{\dot{x} + \dot{y}} = \frac{(\Delta l)^3}{2} \sum\limits_{j = 1}^{n}{\frac{{\partial}z_j}{{\partial}t} \frac{{\partial}\bar{z}_j}{{\partial}t}};z_j = e^{i \theta_j } ; \\ U_{h} = \frac{\Delta l}{2} \sum\limits_{j = 2}^{n} \left(\frac{\theta_j - \theta_{j - 1}}{\Delta l}\right)^2 - F_x \Delta l \sum\limits_{j = 1}^{n}{\cos {\theta_j}} - F_y \Delta l \sum\limits_{j = 1}^{n}{\sin{\theta_j}}. \end{gather*}$$

Здесь появился аналог конечных элементов или схем с центральными разностями. Интегралы в этом выражении заменены конечными изломами, фактически эти интегралы вычислены методом трапеций, т.е. погрешность в определении лагранжиана Lh есть $$O(\Delta l)^{2}$$.

Теперь для построения системы уравнений надо записать

$$$ \frac{d}{dt} \left({\frac{{{\partial}L_{h}}}{{{\partial}\dot \theta_m }}}\right) - \frac{{{\partial}L_{h}}}{{{\partial}\theta_m }} = 0, m = 1, 2, \ldots , n $$$

т. е. продифференцировать дискретный аналог функционала по всем значениям $$\theta_m , \dot \theta_m$$ на введенной сетке. Для рассматриваемой задачи последнее равенство приводит к соотношению

$$$ \sum\limits_{j = 1}^{n}{\frac{{{\partial}^2 L_{h}}} {{{\partial}\theta_k {\partial}\theta_j }} \ddot \theta_j } = \frac{{{\partial}L_{h}}}{{{\partial}\theta_k }} - \frac{{{\partial}^2 L_{h}}}{{{\partial}\dot \theta_k {\partial}t}} - \sum\limits_{j = 1}^{n}{\frac{{{\partial}^2 L_{h}}} {{{\partial}\theta_k {\partial}\theta_j }} \dot \theta_j } . $$$

После подстановки в последнее выражение дискретного аналога лагранжиана и выполнения дифференцирования по всем $$\theta_k , \dot {\theta}_k$$, получаем:

$$\begin{gather*} \sum\limits_{k = 1}^{n}{\alpha_{lk} \ddot {\theta}_k } = n^4 (\theta_{l - 1} - 2 \theta_{l} + \theta_{l + 1} ) + \\ + n^2 (F_y \cos \theta_{l} - F_x {\sin}\theta_{l} ) - \sum\limits_{k = 1}^{n}{g_{lk} {\sin}(\theta_{l} - \theta_k ) \dot {\theta}_k ^2}. \end{gather*}$$

(известно, что $$n\Delta l = 1$$ по построению).

Таким образом, при использовании сеточного аналога вариационного принципа Гамильтона получена дифференциально - разностная система уравнений (дифференциальная по времени, разностная по пространственным переменным).

(рис 9.3)

Несколько затруднительно решать эту систему как систему обыкновенных дифференциальных уравнений, так как она не приведена к нормальной форме Коши. Однако с последней дифференциально - алгебраической системой можно работать и решать ее. Алгоритмы решения основаны на том, что glk — сеточный аналог функции Грина. Обратная матрица — сеточная аппроксимация оператора второй производной. Подробное изложение метода решения в [19.1].

9.2. Вариационные схемы для решения задач газовой динамики

Более интересные разностные схемы получаются для многомерных задач. Особенно плодотворным вариационный подход оказывается при решении задач механики сплошной среды.

Рассмотрим следующий пример. В бесконечном объеме (вакууме) находится область $$\Omega,$$ занятая газом (см. рис. 9.3). Масса газа $$m = \int\limits_\Omega {\rho dV} , \rho (x, y, z)$$ — плотность среды.

Движение системы подчиняется следующим уравнениям:

$$\begin{gather*} \rho \frac{{d{\mathbf{W}}}}{dt} + {grad}{P} = 0 \quad \mbox{(движение)}; \\ \frac{{d \rho }}{dt} + \rho {div}{\mathbf{W}} = 0 \quad \mbox{(неразрывность)}; \\ \rho \frac{{d \varepsilon }}{dt} + P{div}{\mathbf{W}} = 0 \quad \mbox{(закон сохранения энергии)}; \\ F(\rho , P, \varepsilon ) = 0 \quad \mbox{(уравнение состояния)}. \end{gather*}$$

Эти уравнения приведены в лагранжевых переменных. Напомним, что полная производная определена как $$$ \frac{d}{dt} = \frac{{\partial}}{{\partial}t} + ({\mathbf{W}}, \nabla ) $$$. Здесь все обозначения традиционные, через $$\mathbf{W}$$ обозначена скорость движения среды, а через V — элемент объема. Рассматривается случай невязкого нетеплопроводного газа. Запишем лагранжиан для движения среды, представив ее как систему частиц, а затем устремим число частиц в бесконечность.

$$$ L = T - U = \int\limits_\Omega {\frac{{({\mathbf{W}}, {\mathbf{W}})}}{2} \rho dV} - \int\limits_\Omega {\rho \varepsilon dV} = \int\limits_\Omega {\left({\frac{{({\mathbf{W}}, {\mathbf{W}})}}{2} - \varepsilon }\right)dm}. $$$

Движение среды должно доставлять минимум функционалу действия

$$S = \int\limits_0^{t}{L({\tau})d {\tau}}.$$

По закону сохранения массы $$\delta (\rho dV) = \delta (dm) = 0$$. Найдем $$\delta S$$ - вариацию S:

$$\delta S = \int\limits_0^{t}{\delta L({\tau})d {\tau}} = \int\limits_0^{t}{d {\tau}} \int\limits_\Omega {(({\mathbf{W}}, \delta {\mathbf{W}}) - \delta \varepsilon )dm}.$$

Согласно первому началу термодинамики в предположении адиабатичности процесса,

$$$ \delta \varepsilon = \frac{P}{{\rho ^2}} \delta \rho . $$$

Используя кинематическое соотношение

$$$ \delta {\mathbf{W}} = \delta \left({\frac{{d{\mathbf{r}}}}{dt}}\right) = \frac{d}{dt}(\delta {\mathbf{r}}) $$$

получим, что

$$$ \frac{{\delta \rho }}{{\rho ^2}} = \frac{{\delta \left({\frac{1}{V}} \right)}}{{\rho \left({\frac{1}{V}}\right)}} = - \frac{{\delta (dV)}}{{\rho dV}} = - \frac{1}{\rho }{div}(\delta {\mathbf{r}}). $$$

Тогда

$$\begin{gather*} \delta S = - \int\limits_0^{t}{d {\tau}} \int\limits_\Omega {\left({\rho \frac{{d{\mathbf{W}}}}{dt} + {grad}{P}}\right) \delta {\mathbf{r}}dV} + \\ + \left. {\left[{\int\limits_\Omega {({\mathbf{W}}, \delta {\mathbf{r}}) \rho dV}}\right]} \right|_0^{t} + \int\limits_0^{t}{d\tau} \oint\oint\limits_{\partial\Omega }{(\delta {\mathbf{r}}, {\mathbf{n}})Pd {\sigma}} = 0. \end{gather*}$$

Независимые вариации траекторий $$\delta {\mathbf{r}}$$ полагают равными нулю при $$\tau = 0$$ и $$\tau = t$$. Считая также $$\left. {(\delta {\mathbf{r}}, {\mathbf{n}})}\right|_{{\partial}\Omega } = 0$$, получим, что вариация лагранжиана равна нулю на решении системы уравнений Эйлера газовой динамики.

(рис 9.4)

Рассмотрим дискретный аналог функционала действия. Для этого введем в $$\Omega$$ сетку с ячейками, пронумерованными по левому нижнему углу. При отображении $$\Omega$$ на квадрат сетка отображается на равномерную квадратную сетку (см. рис. 9.4). Тогда

$$$ L_{h} = \sum\limits_{i, j}{m_{ij} \left({ \frac{1}{2} < u^2 + v^2 >_{ij} - \varepsilon_{ij}}\right)}, $$$

где $$m_{ij} = \rho _{ij}V_{ij}$$ — масса ячейки ( $$\rho _{ij}$$ — ее плотность) < u2 + v2 > — усредненная по ячейке плотность кинетической энергии

$$$ < u^2 + v^2 >_{ij} = \frac{1}{4} \sum\limits_{l, k = 0}^1 {(u_{i + l, j + k}^2 + v_{i + l, j + k}^2 )}. $$$

Условие адиабатичности течения имеет следствием для дискретной системы выражение $$m_{ij} d\varepsilon _{ij} = - P_{ij} dV_{ij}$$.

Кинематические соотношения (уравнения движения для узлов сетки) есть

$$$ \frac{dx_{ij}}{dt} = u_{ij}, \frac{dy_{ij}}{dt} = v_{ij}. $$$

Объем каждой ячейки вычисляется в предположении, что границы ячейки — отрезки прямых:

$$$ V_{ij} = \frac{1}{2}[(x_{i + 1j} - x_{ij + 1})(y_{i + 1j + 1} - y_{ij} ) - (x_{i + 1j + 1} - x_{ij} )(y_{i + 1j} - y_{ij + 1} )]. $$$

Теперь необходимо получить соответствующую систему дифференциально - разностных соотношений, аппроксимирующих уравнения Эйлера.

Запишем функционал действия

$$S_{h} = \int\limits_0^{t}{L_{h} ({\tau})d {\tau}}$$

и найдем его вариацию:

$$\begin{gather*} 0 = \delta S_{h} = \int\limits_0^{t}{\delta \left({\sum\limits_{\omega_{h}}{m_{ij} \left({\frac{{< u^2 + v^2 >_{ij}}}{2} - \varepsilon_{ij}}\right)}}\right)} d {\tau} = \\ = \sum\limits_{\omega_{h}}{\left\{{\delta \left({\frac{{m_{ij}}}{8} \sum\limits_{l, k = 0}^1 {(u_{i + l, j + k}^2 + v_{i + l, j + k}^2 )}}\right) - \delta (m_{ij} \varepsilon_{ij} )}\right\}} = \\ = \sum\limits_{\omega_{h}}{\frac{1}{8} \delta (m_{ij} \sum\limits_{l, k = 0}^1 {u_{i + l, j + k}^2} )} + \sum\limits_{\omega_{h}}{\frac{1}{8} \delta (m_{ij} \sum\limits_{l, k = 0}^1 {v_{i + l, j + k}^2} )} + \sum\limits_{\omega_{h}}{P_{ij} \delta V_{ij}}. \end{gather*}$$

Первое слагаемое в последнем равенстве можно переписать как

$$$ \sum\limits_{\omega_{h}}{\delta (u_{ij}^2 \cdot \frac{1}{8} \sum\limits_{l, k = 0}^1 {m_{i - k, j - l}} )} = M_{ij} u_{ij} \delta u_{ij} = M_{ij} \frac{{du_{ij}}}{dt} \delta x_{ij}, $$$

где

$$$ M_{ij} = \frac{1}{4}(m_{i - 1j} + m_{i - 1j - 1} + m_{ij} + m_{ij - 1} ). $$$

После преобразований получаются следующие соотношения:

$$\begin{gather*} M_{ij} \frac{{du_{ij}}}{dt} = \sum\limits_{l, k = 0}^1 {P_{i - l, j - k} \frac{{{\partial}V_{i - l, j - k}}}{{{\partial}x_{ij}}}} ; \\ M_{ij} \frac{{dv_{ij}}}{dt} = \sum\limits_{l, k = 0}^1 {P_{i - l, j - k} \frac{{{\partial}V_{i - l, j - k}}}{{\partial y_{ij}}}} ; \end{gather*}$$

т.е. возникает система дифференциально - разностных уравнений движения.

Уравнение адиабатичности трактуется, как уравнение для определения энергии:

$$\begin{multline*} M_{ij} \frac{{d \varepsilon_{ij}}}{dt} = - P_{ij} \frac{{dV_{ij}}}{dt} = \\ = - P_{ij} \sum\limits_{l, k = 0}^1 {\left({\frac{{\partial V_{ij}}}{{{\partial}x_{i + l, j + k}}} u_{i + l, j + k} + \frac{{{\partial}V_{ij}}}{{{\partial}y_{i + l, j + k}}}y_{i + l, j + k}}\right)} ; \end{multline*} $$$

Система дифференциально - разностных соотношений замыкается уравнением состояния.

Теперь, заменяя производные по времени на конечные разности, получим разностную схему для решения уравнений газовой динамики на нерегулярной подвижной сетке. Обычно для решения таких уравнений используются явные разностные схемы. Подробнее о вариационных методах получения схем для уравнений газовой динамики в [19.2].

9.3. Вариационная схема для уравнения теплопроводности на криволинейной сетке

(рис 9.5)

Рассмотрим линейное уравнение теплопроводности

$$$ {\frac{{\partial}u}{{\partial}t} + {div}(- k \ {grad}\ {u}) = 0} $$$

с условиями

$$(-k \ {grad}\ {u}, {\mathbf{n}}) = 0$$

в ограниченной области $$\Omega$$ с криволинейной границей. При этом нигде в $$$ \bar {\Omega} $$$ уравнение не вырождено, т.е. k(x, y) > 0 во всех точках области, включая граничные.

В области $$\Omega$$ каким - либо образом введена сетка с четырехугольными ячейками. Сетка считается связной, т.е. для любых двух вершин ячеек существует ломаная, их соединяющая и состоящая из ребер ячеек (рис. 9.5).

Пусть сетка построена так, что существует преобразование, переводящее область $$\Omega$$ в параллелограмм (прямоугольник) с равномерной сеткой внутри. Тогда координатные линии x, y переходят в координатные кривые криволинейного базиса $$\xi , \eta$$.

Перепишем уравнение (9.1) в виде системы

$$\begin{gather*} \frac{{\partial}u}{{\partial}t} + div{\mathbf{W}} = 0, \\ {\mathbf{W}} + k {grad} u = 0. \end{gather*}$$

Рассмотрим функционал

$$$ {F[u] = \int\limits_\Omega {\left({\frac{{({\mathbf{W}}, {\mathbf{W}})}}{k} - \frac{{\partial}}{{\partial}t}u^2}\right)dxdy}.} $$$

Найдем $$\delta F[u]$$:

$$\begin{gather*} \delta F[u] = \int\limits_\Omega {\delta \left({\frac{{({\mathbf{W}}, {\mathbf{W}})}}{k} - \frac{{\partial}}{{\partial}t}u^2}\right)dxdy} = \\ = \int\limits_\Omega {\left({- 2{div}{\mathbf{W}} \cdot \delta u - 2 \delta u \frac{{\partial}u}{{\partial}t} - 2u \frac{{\partial}}{{\partial t}} \delta u}\right)dxdy}, \\ \delta \frac{{({\mathbf{W}}, {\mathbf{W}})}}{k} = 2 \left({\frac{{\mathbf{W}}}{k}, \delta {\mathbf{W}}}\right) = 2 \left({\frac{{\mathbf{W}}}{k}, - k {grad} \delta u}\right) = - 2(\delta u, {div}{\mathbf{W}}). \end{gather*}$$

Отсюда при $$$ \frac{{\partial}}{{\partial}t} \delta u = 0 $$$ минимум функционала достигается на решении уравнения теплопроводности.

Для построения разностной схемы введем дискретный аналог функционала $$F_{h}({\mathbf{W}}^{h} )$$, т.е. в дискретном аналоге основной расчетной величиной будет поток тепла.

Прежде чем построить функционал, рассмотрим ячейку разностной сетки (рис. 9.6). Температуру uij и коэффициент теплопроводности (или температуропроводности) kij отнесем к центру ячейки (точке пересечения диагоналей). В дальнейшем считаем, что термодинамические величины постоянны во всей ячейке. Векторы теплового потока отнесем к углам ячейки (рис. 9.6), а к центрам соответствующих ребер — проекции потоков на координатные оси. Считаем, что i увеличивается по мере увеличения координаты $$\xi$$ ; j — по мере увеличения $$\eta$$ ; проекции векторов потока направлены вдоль соответствующих координатных линий. Заметим, что проекции потоков для двух ячеек сонаправлены с векторами внешней нормали, а для двух — противонаправлены.

(рис 9.6)

Проинтегрируем уравнение

$$$ \frac{{\partial}u}{{\partial}t} + {div}{\mathbf{W}} = 0 $$$

по элементарной ячейке разностной сетки. Имеем:

$$$ {S_{ij} \frac{{du_{ij}}}{dt} + (W_{\eta ij} \Gamma_{\xi ij} - W_{\xi i + 1j} \Gamma_{\eta i + 1j} - W_{\eta ij + 1} \Gamma_{\xi ij + 1} + W_{\xi ij} \Gamma_{\eta ij} ) = 0, } $$$

где $$\Gamma$$ — длины соответствующих ребер, Sij — площадь элементарной ячейки.

Так как координаты всех вершин выпуклого четырехугольника известны, то поиск длин, площадей и углов — элементарная геометрическая задача.

Уравнение (9.4}) — дискретный аналог уравнения (9.2). Если возможно определить все потоки в моменты времени tn, tn + 1, а после применить аппроксимацию (9.2) по времени с какими - либо весами, то будет построена разностная схема для расчета температуры.

Учтем, что

$$$ 2u \frac{{\partial}u}{{\partial}t} = - 2u{div}{\mathbf{W}}. $$$

Построим дискретный аналог (9.3):

$$\begin{gather*} F_{h} ({\mathbf{W}}^{h} ) = \sum\limits_{i, j \in \omega_{h}} \left[{S_{ij} \cdot \left(\sum\limits_{p, l = 0}^1{\frac{{W_{i + p, j + l}^2}}{{k_{ij}}}}\right) +}\right. \\ {\left. \begin{array}{l} {} \\ {} \\ {} \\ \end{array}2u_{ij}^{({\sigma})} (W_{\eta ij} \Gamma_{\xi ij} - W_{\xi i + 1j} \Gamma_{\eta i + 1j} - W_{\eta ij + 1} \Gamma_{\xi ij + 1} + W_{\xi ij} \Gamma_{\eta ij}) \right].}\end{gather*}$$

Скалярные квадраты, входящие в первое слагаемое дискретного аналога функционала, выражаются через контравариантные проекции следующим образом (рис. 9.7):

$$\begin{gather*} W_{ij} = \frac{1}{{{\sin}^2 \varphi_1}}(W_{\eta ij}^2 + W_{\xi ij}^2 + 2W_{\eta ij} W_{\xi ij} \cos \varphi_1 ), \\ W_{i + 1j} = \frac{1}{{{\sin}^2 \varphi_2}}(W_{\eta ij}^2 + W_{\xi i + 1j}^2 - 2W_{\eta ij} W_{\xi i + 1j} \cos \varphi_2 ) \end{gather*}$$

и т.д.

(рис 9.7)

Знак " + " или " - " определяется по правилу: компоненту потока приписывается знак " + ", если проекция потока сонаправлена с внешней нормалью, а знак " - " — если противонаправлена. Таким образом, для углов ячейки $$\varphi_1$$ и $$\varphi_3$$ получим знак " + " в последнем слагаемом (проекции одинаковых знаков), а для углов $$\varphi_2$$ и $$\varphi_4$$ — знак " - " (проекции теплового потока в произведении разных знаков).

Для получения явной схемы положим в (9.5) вес верхнего слоя по времени $$\sigma = 0$$ и дифференцируем (9.5) по всем $$W_{\xi ij}, W_{\eta ij}$$. Приравнивая производные нулю, получим схему для определения потоков, затем из (9.4) ищем все $$u_{ij}^{{n} + 1} $$.

Для построения неявной схемы в (9.4) считаем $$\sigma = 1$$, а вместо (9.5) пишем следующую дискретизацию:

$$$ S_{ij} \frac{{u_{ij}^{{n} + 1} - u_{ij}^{n}}}{\tau} = - \sum\limits_{l, p = 0}^1 {(- 1)^{l}(W_\xi )_{i + l, j}^{{n} + 1} \Gamma_{\eta i + l, j} + (- 1)^{p} (W_\eta )_{i, j + p}^{{n} + 1} \Gamma_{\xi i, j + p} .} $$$

Выражая отсюда неизвестное пока значение $$u_{ij}^{{n} + 1}$$ в (9.5), получим выражение, зависящее от $$u_{ij}^{n}, \{W_\xi ^{{n} + 1} \}, \{W_\eta ^{{n} + 1} \}$$, причем $$F_{h} ({\mathbf{W}}^{h} )$$ есть сумма квадратов контравариантных проекций.

Дифференцируя, получим линейную систему уравнений для определения потоков. Можно показать, что матрица системы будет обладать следующими свойствами:

$${\mathbf{A}} = {\mathbf{A}}^* > 0;$$

$$\mathbf{A}$$ имеет ленточную структуру;

$$\mathbf{A}$$ является разреженной.

Можно применить эффективные итерационные методы решения системы.

Доказано, что неявная схема будет безусловно устойчивой, а явная — условно устойчивой.

Метод легко обобщается на случай k = k(x, y, u), если уравнение не вырождается. Кроме того, метод может быть обобщен и на случай других граничных условий (не обязательно отсутствия потоков). В этом случае в функционал (9.3) добавляются соответствующие интегралы по границам, а в (9.5) — суммы по поверхностям.

Подробнее об этих схемах можно прочитать в [19.2].

9.4. Задачи для самостоятельного решения

Уравнение Кортевега - Де Фриза

Одно из самых замечательных уравнений математической физики — уравнение Кортевега - Де Фриза (сокращенно КДФ) часто записывают в виде

ut - 6uux + uxxx = 0

или

$$$ \tilde {u}_t + \tilde {u} \tilde {u}_x = \tilde {u}_{xxx}. $$$

Найти преобразование, переводящее эти формы записи друг в друга.

Рассматриваем задачу для уравнения ut - 6uux + uxxx = 0 в области $$x \in [- 10;10]$$ с условием периодичности.

Для решения иногда используют трехслойную разностную схему на шаблоне рис. 9.8 (третья производная расписывается по пяти точкам симметричным образом, некоторые коэффициенты могут обратиться в нуль). Исследовать ее на аппроксимацию и устойчивость. Какое условие устойчивости получено? Построить разностную схему на шаблоне, напоминающем шаблон схемы Саульева для решения уравнения теплопроводности (рис. 9.9 а, б).

(рис 9.9) (рис 9.8)

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

Уравнение ut - 6uux + uxxx = 0 имеет бесконечное число законов сохранения. Укажем несколько из них:

$$\begin{gather*} \int {udx = {const}_{1} , } \\ \int {u^2 dx = {const}_{2}} \\ \int {\left({\frac{{(u^{\prime}_{x})^2}}{2} + u^3 }\right)} dx = {const}_{3} = I_1 . \end{gather*}$$

Как построить консервативную разностную схему, чтобы на сеточном уровне выполнялись законы сохранения $$\int {udx} = {const}{?}$$

Особую роль играет третий из приведенных выше законов сохранения. Он является гамильтонианом для уравнения ut - 6uux + uxxx = 0, т.е.

$$$ \frac{{\partial}u}{{\partial}t} = \frac{d}{dx} \left({\frac{{\delta I_1}}{{\delta u}}}\right). $$$

Здесь $$$ \frac{\delta I_1}{\delta u} $$$ вариационная производная функционала I1, способ получения разностной схемы, сохраняющей гамильтониан системы.

Решение. Запишем сеточный аналог гамильтониана I1:

$$$ I_{{1}}^{{h}} = \sum\limits_{m = - \infty }^{+ \infty } {\left[{\frac{1}{2} \left({\frac{{u_{{{m}} + {{1}}} - u_{{m}}}}{h}}\right)^2 + (u_{{m}} )^3 }\right]h = I_1 + O(h^2 )}. $$$

Контрольный вопрос: почему $$I_1^{h}$$ аппроксимирует I1 с точностью O(h2)?

На сеточном уровне взятие вариационной производной означает дифференцирование по всем $$u_{m}^{n}$$ и деление на h. Тогда получаем сеточную запись вариационной производной:

$$\begin{gather*} \frac{\delta I_1^h}{\delta u} = \frac{1}{h}\frac{{\partial}_1^h}{{\partial}u_m } = u_m^2 + \frac{1}{h^2}((u_m - u_{m - 1}) - (u_{m + 1} - u_m )) = \\ = 3u_m^2 - \frac{u_{m + 1} - 2u_{m} + u_{m - 1}}{h^2} + O(h^2). \end{gather*}$$

Аппроксимируя дискретный аналог (4.4.3) с естественным для этого вторым порядком по h, получаем

$$\frac{{\partial} u_m}{{\partial} t} = \frac{3(u_{m + 1})^2 - 3(u_{m - 1})^2}{2h} - \frac{u_{m + 2} - 2u_{m + 1} + 2u_{m - 1} - u_{m - 2}}{2h^3} .$$

Заменяя производную по времени разностью$$\frac{u_m^{n + 1} - u_m^{n - 1}}{2 {\tau}}$$ и вычисляя правую часть на n слое, получаем одну из схем пункта 2.

Конечно, возможны и другие аппроксимации гамильтониана, варьирование которых приводит к другим разностным схемам. Все они будут записываться на симметричных шаблонах, на сеточном уровне для этих схем также будет выполняться закон сохранения $$\int {udx = {const}}$$.

Удастся ли получить вариационную схему на несимметричном шаблоне типа Саульева?