Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 1279 / 290 | Оценка: 4.40 / 4.36 | Длительность: 21:57:00
Специальности: Математик
Лекция 3:

Численные методы решения уравнений в частных производных гиперболического типа (на примере уравнения переноса)

  1. Исследовать на спектральную устойчивость разностную схему на шаблоне квадрат

    \begin{gather*}  \frac{1}{{2{\tau}}} \left[{(u_{m + 1}^{n + 1} + u_m^{n + 1}) - (u_{m + 1}^{n} - u_m^{n})}\right] + \\ 
 + \frac{1}{h} \left[{(u_{m + 1}^{n + 1} + u_{m + 1}^{n}) - (u_m^{n + 1} - u_m^{n})}\right] = f_{m + 1/2}^{n + 1/2},    \end{gather*}

    аппроксимирующую линейное уравнение переноса

    $ \frac{{\partial}u}{{\partial}t} + \frac{{\partial}u}
{{\partial}x} = 0.  $

    Решение. После подстановки в схему гармоник Фурье

    u_m^{n} = {\lambda}^{n} e^{i{\alpha} m},

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

    $ {\lambda}(\alpha ) = \frac{1 - i{\sigma}{\tg{\alpha}/2}}{1 + i{\sigma}{\tg{\alpha}/2}}, {\sigma}={\tau}/h \mbox{ - число Куранта.}  $

    Тогда \left| {{\lambda}(\alpha )}\right| = 1 и рассматриваемая схема безусловно устойчива. Этого следовало ожидать, так как аппроксимация производной по времени производится как на нижнем, так и на верхнем временном слое. Однако решение на верхнем слое выписывается в явном виде, поскольку правое условие ставится на левой границе. Поэтому схему иногда относят к явным схемам бегущего счета. Расчетная формула будет

    $ u_{m + 1}^{n + 1} = u_m^{n} +  \frac{{{\sigma}- 1}}{{{\sigma}+ 1}}(u_m^{n + 1} - u_{m + 1}^{n}) + {\tau}f_{m + 1/2}^{n + 1/2} .  $

    Схема обладает вторым порядком аппроксимации по времени и пространственной координате.

  2. Получить дисперсионное соотношение для разностных схем первого и второго порядка сходимости для аппроксимации нелинейного уравнения

    $ \frac{{\partial}v}{{\partial}t} + u \frac{{\partial}v}{{\partial}x} = 0, \quad \mbox{где} \quad u \ne {const},  $

    и, вообще говоря, зависит от решения.

    Указание. Использовать дисперсионное соотношение для данного дифференциального уравнения в виде e^{\lambda t} e^{ikx}, где k — волновое число, \lambda  = - uki.

    Решение. Будем искать дисперсионное соотношение для разностного уравнения в виде

    v_m^{n} = e^{{\lambda}(k)t_n } \cdot e^{{ikx}_m } = e^{{\lambda}n{\tau}+ {ikmh}}.

    После его подстановки в схему правый уголок первого порядка аппроксимации, получим

    $ \frac{1}{\tau} \left({e^{{\lambda}{\tau}} - 1}\right) + 
 \frac{u}{h} \left({1 - e^{- {ikh}}}\right) = 0,  $

    откуда следует

    $ {\lambda}\left({{\tau}, h, k}\right) ={\tau}^{- 1} \ln \left\{{1 - 
 \frac{u \tau}{h} + \frac{u \tau}{h}e^{- ikh}}\right\}.  $

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

    Пусть

    $ \frac{{u{\tau}}}{h} = 1.  $
    Имеет место совпадение \lambda (k) и \lambda (\tau , h, k), однако этот случай мало интересен для вычислительной практики, поскольку u может меняться в ходе решения задачи. Для проведения дисперсионного анализа схем положим kh \ll 1, имея в виду то, что h — малый параметр, а аппроксимация тем лучше, чем меньше волновое число, т.е. чем более гладким является исследуемое частное решение. Заметим, что для расчетной сетки с шагом h обычно реализуются волновые числа k  \approx  2 \pi h^{- 1}.

    В таком случае дисперсионное соотношение будет

    \begin{gather*}  {\lambda}\left({{\tau}, {h, k}}\right)  \approx  - {iuk} - 
 \frac{{{uhk}^2}}{2} \left({1 - \frac{{u{\tau}}}{h}}\right) = {\lambda}(k) - \frac{{k^2}}{2}{uh} \left({1 - {\sigma}}\right), \\ 
 {\sigma}= u{\tau}/h \mbox{ - число Куранта, }   \end{gather*}

    а частное решение примет вид

    $ \exp (ik({mh} - {un}{\tau})) \cdot \exp \left({-  \frac{1}{2}{uk}^2 (h - 
u{\tau})n{\tau}}\right).  $

    Первый сомножитель совпадает с соответствующим частным решением дифференциального уравнения, второй является особенностью разностного решения; его поведение при различных \tau , h, k представляет особенный интерес. При u < 0, h - u\tau  > 0, как видно из этого решения, рассматриваемую схему нельзя использовать для проведения расчетов. Второй сомножитель имеет порядок exp (k2 | u | htn), он быстро растет при k \sim h^{- 1}, h \to 0,{\tau}\to 0. Аналогичная ситуация получается и при u > 0, h - u\tau  < 0. При u < 0, h - u\tau  > 0 второй сомножитель затухает с ростом tn тем быстрее, чем больше k, или чем меньше \lambda — длина волны частного решения вида exp(ikx). Таким образом, численное решение, полученное по схеме уголок, отличается от точного, в котором все гармоники сохраняют амплитуду с ростом tn, наличием затухающего множителя для гармоник с большими волновыми числами (или малыми длинами волн). Действие этого сомножителя приводит к сглаживанию решений, имеющих разрыв в начальных данных.

    Для разностной схемы Лакса - Вендроффа второго порядка

    $ (v_{m + 1}^{n + 1} + v_m^{n}) - \sigma (v_m^{n} - v_{m - 1}^{n}) + \frac{{\sigma}}{2} (v_{m - 1}^{n} - 2v_m^{n} + v_{m + 1}^{n}) = 0  $

    дисперсионное соотношение будет

    $  {\lambda} ={\tau}^{- 1} \ln \left\{{1 + {\sigma}(e^{- {ikh}} - 1) + {2uh}(1 -   \sigma) {\rm sin}}^2 \frac{{kh}}{2}\right\},    $

    для длинноволновых гармоник с kh \ll {1}

    $ {\lambda}\left({{\tau}, {h, k}}\right)  \approx  - {iku} + i \frac{{k^3 u}}{6} (h^2 - u^2{\tau}^2 ) = {\lambda}(k) + ik^3 \frac{{{uh}^2}}{6} (1 - {\sigma}^2 ).  $

    Из полученных соотношений можно сделать следующий вывод. Решения исходного уравнения переноса имеют вид волн с волновым числом k, которые движутся вправо со скоростью u, так как

    v(t, x) = exp(ikx + \lambda (k)t) =  exp(ik(x - ut)).

    Решения разностного уравнения имеют вид

    \begin{gather*}  v \left({t_n , x_m }\right) = \exp ({ikx}_m + 
 {\lambda}({\tau}, {h, k})t_n ) = \\ 
 = \exp \left({ik \left\{{x_m - u \left[{1 - \frac{{k^2}}{6} \left({h^2 - u^2{\tau}^2}\right)}
\right]t_n }\right\}}\right) = \\ 
 = \exp \left({ik \left\{{x_m - u \left[{1 + {U}_k (k)}\right]t_n }\right\}}\right),  
\end{gather*}

    где введено обозначение

    $ U_k = - \frac{k^2}{6} \left({h^2 - u^2{\tau}^2}\right).  $

    Последнее выражение означает, что разностные решения уже не имеют вид волн, движущихся с одной и той же скоростью. Теперь каждая волна со своей частотой движется с собственной скоростью u = Uk[1 + u(k)]. Разумеется, при малых k и скорости uk мало отличаются от u, но высокочастотные волны движутся со скоростями, заметно отличающимися от скорости переноса. Кроме того, в схеме Лакса - Вендроффа гармоники со временем не затухают, численное решение не сглаживается со временем. Начальный волновой пакет u и определенный им начальный профиль решения, изменяются с течением времени из - за рассогласования фаз. Это приводит к потере монотонности профиля v(x), если он был вначале монотонным, и появлению осцилляций разностного происхождения. Появление сеточной дисперсии — одно из проявлений эффекта Гиббса, подробнее об эффекте Гиббса в [13.16].

  3. Акустическая система описывает распространение плоских звуковых волн. Поставим для нее задачу Коши:
    \begin{gather*}  \frac{{\partial}u}{{\partial}t} + \frac{1}{\rho} \frac{{\partial}p}{{\partial}x} = 0, \frac{{\partial}p}{{\partial}t} + {\rho}c^2 \frac{{\partial}u}{{\partial}x} = 0, \\ 
 {\rho} = {const}, \quad c = {const}, \\ 
u(x, 0) = u_0 (x), \quad p(x, 0) = p_0 (x), \\ 
t > 0, - \infty  < x < \infty ,  \end{gather*}

    где u — скорость движения среды, p — давление, \rho - плотность среды, c — скорость звука в среде.

    • a) Представить акустическую систему в интегральной форме.
    • b) Преобразовать акустическую систему к виду системы с разделенными переменными.
    • c) Предложить разностную схему первого порядка точности для ее решения.

    Решение.

    • a) Проинтегрировав акустическую систему по произвольной области с границей G, получим:

      \begin{gather*}  \oint\limits_{\Gamma}{{\rho}udx - {pdt}} = 0, \\ 
 \oint\limits_{\Gamma}{\frac{p}{c^2}dx - {\rho}{udt}} = 0.  \end{gather*}

    • b) Умножим второе уравнение на (\rho c)^{ - 1}, затем сложим с первым и вычтем из него. Получим систему двух уравнений

      \begin{gather*}
 \frac{{\partial}}{{\partial}t} \left({u + \frac{p}{{\rho}c}}\right) + c \frac{{\partial}}{{\partial}x} \left({u + \frac{p}{{\rho}c}}\right) = 0 \\ 
 \frac{{\partial}}{{\partial}t} \left({u - \frac{p}{{\rho}c}}\right) - c \frac{{\partial}}{{\partial}x} \left({u - \frac{p}{{\rho}c}}\right) = 0, 
 \end{gather*}
  $

      или, после введения обозначений (инвариантов Римана)

      $ R = u + \frac{p}{{\rho}c}, S = u - \frac{p}{{\rho}c}  $
      система запишется в виде двух уравнений переноса:

      $ \frac{{{\partial}R}}{{\partial}t} + c \frac{{{\partial}R}}{{\partial}x} = 0 \quad \frac{{{\partial}S}}{{\partial}t} - c \frac{{{\partial}S}}{{\partial}x} = 0.  $

    Эта система позволяет выписать общее решение:

    R = f(x - ct),  S = g(x + ct).

    Здесь f, g — произвольные непрерывно дифференцируемые функции, определяемые из начальных условий. В терминах инвариантов Римана можно определить и значения естественных переменных u, p:

    \begin{gather*}  u = \frac{1}{2}(f(x - ct) + g(x + ct)), \\ 
p = \frac{{\rho}c}{2}(f(x - ct) - g(x + ct)).  \end{gather*}

    Замечательное свойство инвариантов Римана R, S заключается в их постоянстве вдоль прямых x - ct = const, и x + ct = const или

    $ \frac{dx}{dt} =  \pm  c  $
    характеристик акустической системы.

    Если задавать начальные данные u(x, 0) = {\varphi}(x) p(x, 0) = {\psi}(x), {- \infty  < x < \infty}, то связь начальных данных и инвариантов Римана будет

    $ {\varphi}\left(x\right) = \frac{1}{2}(f(x) + g(x)),  {\psi}(x) = \frac{{\rho}c}{2} \left[{f(x) - g(x)}\right],  $

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

    $ f(x) = {\varphi}(x) + \frac{{{\psi}(x)}}{{\rho}c}, g(x) = {\varphi}(x) - \frac{{{\psi}(x)}}{{\rho c}}.  $

    Следовательно, решение системы в этом случае

    \begin{gather*}
u(t, x) = \frac{1}{2}({\varphi}(x - ct) + {\varphi}(x + ct)) +  \frac{1}{{2 {\rho}c}}({\psi}(x - ct) + {\psi}(x + ct)), \\ 
p(t, x) = \frac{1}{2}({\varphi}(x - ct) + {\varphi}(x + ct)) -  \frac{{\rho}c}{2}({\psi}(x - ct) + {\psi}(x + ct)).
 \end{gather*}

    Для численного решения системы, записанной в инвариантах Римана, можно применить схему Куранта - Изаксона - Риса:

    \begin{gather*}  \frac{{R_m^{n + 1} - R_m^{n}}}{\tau} + 
c \frac{{R_m^{n + 1} - R_{m - 1}^{n}}}{h} = 0, \\ 
 \frac{{S_m^{n + 1} - S_m^{n}}}{\tau} - c \frac{{S_{m + 1}^{n} - S_m^{n}}}{h} = 0, 
 \end{gather*}

    или

    \begin{gather*}  
R_m^{n + 1} = \left({1 - {\sigma}}\right)R_m^{n} + {\sigma}R_{m - 1}^{n},  \\ 
S_m^{n + 1} = \left({1 + {\sigma}}\right)S_m^{n} - {\sigma}S_{m + 1}^{n} .  
\end{gather*}

    Известно, что данная схема имеет первый порядок аппроксимации и является устойчивой при выполнении условия КФЛ:

    $ \frac{{c{\tau}}}{h} \le 1.  $

  4. Заменить смешанную задачу для волнового уравнения

    \begin{gather*}  \frac{{{\partial}^2 u}}{{\partial}t} = c^2 \frac{{{\partial}^2 u}}{{{\partial}x^2}} + f(t, x), \quad 0  <  t {\le} T, \quad 0  <  x  <  X, u(0, x) = \varphi_1 (x), \\ 
 \frac{{{\partial}u(0, x)}}{{\partial}t} = \varphi_2 (x), 0  <  x  <  X, u(t, 0) = \varphi_3 (t), 
u(t, X) = \varphi_4 (t), 0 \le t \le T   \end{gather*}

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

    Решение. Введем новые переменные

    $ v(t, x) = \int\limits_0^{x}{\frac{{{\partial}u(t, \eta )}}{{\partial}t}d \eta }, F(t, x) = \int\limits_0^{x}{f(t, \eta )d \eta } .  $

    Функции u(t, x), v(t, x), что несложно проверить, удовлетворяют акустической системе уравнений

    \begin{gather*}  \frac{{\partial}u}{{\partial}t} - \frac{{\partial}v}
{{\partial}x} = 0, \frac{{\partial}v}{{\partial}t} - c^2 \frac{{\partial}u}{{\partial}x} = F(t, x), \\ 
u(0, x) = \varphi_1 (x), v(0, x) = \int\limits_0^{x}{\varphi_2} (\eta )d \eta , u(t, 0) = \varphi_3 (t),  u({t, X}) = \varphi_4 (t).  \end{gather*}

    Введем в рассмотрение параметрическую схему для численного решения этой системы:

    \begin{gather*}
 \frac{{u_m^{n + 1} - u_m^{n}}}{\tau} =  \frac{{c^2}}{h} \left\{{\xi_1 (v_{m + 1}^{n + 1} - 
v_m^{n + 1}) + (1 - \xi_1 )(v_{m + 1}^{n} - v_m^{n})}\right\} + F_{m + 1/2}^{n + 1/2}, \\ 
 \frac{{v_m^{n + 1} - v_m^{n}}}{\tau} =  \frac{1}{h} \left\{{\xi_2 (u_m^{n + 1} - u_{m - 1}^{n + 1}) + (1 - \xi_2 )(u_m^{n} - u_{m - 1}^{n})}\right\}, \\ 
0 \le \xi_1 , \xi_2 \le 1.
 \end{gather*}

    Исследование полученной двухпараметрической разностной схемы на спектральную устойчивость дает условия устойчивости:

    $ \xi_1 + \xi_2 \ge 1, {\sigma}^2 \left({2 \xi_1 - 1}\right) \left({2 \xi_2 - 1}\right) \ge - 1, {\sigma}= \frac{{c \tau }}{h}.  $
    Отсюда видно, что если \xi_1 \ge 0, 5 и {\xi}_2 \ge 
0, 5, то схема устойчива при любых числах Куранта и сеточных параметрах, если \xi_1 + {\xi}_2 \ge 1, но один из параметров меньше 0, 5, то схема условно устойчива при
    ${\tau}\le \frac{h}{c}(2 \xi_1 - 1)^{- 1} (2 \xi_2 - 1)^{- 1}.  $

    Пусть \xi _{1} = \xi _{2} = 0, 5. Тогда полученная схема является безусловно устойчивой и обладает вторым порядком аппроксимации по обоим переменным.

  5. Исследовать на сходимость схему "крест" с шаблоном

    для численного решения волнового уравнения

    $ \frac{{\partial}^2 u}{{\partial}t^2} -  \frac{{\partial}^2 u}{{\partial}x^2} = 0.  $

    Решение. Запишем разностную схему:

    $ {\mathbf{L}}_{\tau} u_{\tau} = \frac{{u_m^{n + 1} - 2u_m^{n} + 
u_m^{n - 1}}}{{{\tau}^2}} - \frac{{u_{m - 1}^{n} - 2u_m^{n} + {u}_{m + 1}^{n}}}{{h^2}} = 0.   $

    Исследование на аппроксимацию дает выражение для главного члена невязки

    $ r_{\tau} = \frac{{{\tau}^2}}{{12}}u_t^{(4)} - \frac{{h^2}}{{12}}
u_x^{(4)} + O({\tau}^2 , h^2 ).  $

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

    $ {\lambda}^2 - 2 \left({1 - 2 {\sigma}^2 {\sin}^2 \frac{\alpha }{2}}
\right) {\lambda}+ 1 = 0,  {\sigma}={\tau}/h \mbox{ - число Куранта.}  $

    По теореме Виета, произведение корней \lambda _{1} \lambda _{2} = 1, и условие устойчивости выполняется, если | \lambda _{1} | = | \lambda _{2} | = 1. Для полученного квадратного уравнения с действительными коэффициентами это означает, что его корни являются комплексно сопряженными. Это возможно лишь в случае, если дискриминант

    $  D(\alpha ) = 4 {\sigma}^2 {\sin}^2{\alpha} \left({{\sigma}^2 {\sin}^2 
\frac{\alpha }{2} - 1}\right)  $

    будет отрицательным. Условие устойчивости будет выполнено для всех гармоник, если \sigma  < 1.