Дискуссионный математический форумМатематический форум
Математический форум Math Help Planet

Обсуждение и решение задач по математике, физике, химии, экономике

Теоретический раздел
Часовой пояс: UTC + 3 часа [ Летнее время ]
новый онлайн-сервис
число, сумма и дата прописью

Часовой пояс: UTC + 3 часа [ Летнее время ]


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

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


Метод прямых


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


Основную идею конструирования метода прямых рассмотрим на примере нахождениия решения u(x,y) уравнения:


\begin{gathered}A(x,y)\frac{\partial^2 u(x,y)}{\partial x^2}+ 2B(x,y)\frac{\partial^2 u(x,y)}{\partial y\,\partial  x}+ C(x,y) \frac{\partial^2 u(x,y)}{\partial y^2}+ D(x,y)\frac{\partial u(x,y)}{\partial x}\,+\\[2pt] +\,E(x,y)\frac{\partial u(x,y)}{\partial y}+ F(x,y)u(x,y)= G(x,y),\quad (x,y)\in \Omega= (0,L)\times (0,M), \end{gathered}
(8.59)

с краевыми условиями


\begin{aligned}& (x,0)= \varphi(x), & & u(x,M)= \psi(x),\\ & (0,y)= \Phi(y), & & u(L,y)= \Psi(y) \end{aligned}
(8.60)

на границах области прямоугольной формы (рис. 8.15), где \varphi(x),\psi(x), \Phi(x), \Psi(x) — заданные функции. Пусть выполнено условие B^2-AC<0 эллиптичности уравнения для всех (x,y)\in \Omega. Отсюда следует, что A(x,y)\ne 0,~ C(x,y)\ne 0.


Нахождение решения ДУ в частных производных на границах области

Для решения задачи выберем на отрезке [0,M] точки y_{j}= j\,h,~ 0 \leqslant j \leqslant J,~ J\,h=M, и проведем через них прямые, параллельные оси Ox. На каждой прямой дифференциальное уравнение в частных производных (8.59) приближенно заменим обыкновенным дифференциальным уравнением для функций u(x,y_{j}). Для этого аппроксимируем частные производные по y\colon


\begin{gathered}\left.{\frac{\partial \widehat{u}(x,y)}{\partial y}}\right|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-\widehat{u}(x,y_{j-1})}{2h}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y^2}}\right|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-2\,\widehat{u}(x,y_{j})+ \widehat{u}(x,y_{j-1})}{h^2}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y\, \partial x}}\right|_{y=y_{j}}= \frac{\widehat{u}_{x}(x,y_{j+1})-\widehat{u}_{x}(x,y_{j-1})}{2h}\,. \end{gathered}
(8.61)

Введем обозначения


\begin{gathered}\widehat{u}(x,y_{j})= \widehat{u}_{j}(x),\qquad \frac{\partial \widehat{u}(x,y_{j})}{\partial x}= \widehat{u}\,'_{j}(x),\qquad \frac{\partial^2 \widehat{u}(x,y_{j}) }{\partial x^2}= \widehat{u}\,''_{j}(x),\\[4pt] A(x,y_{j})= A_{j}(x),\qquad B(x,y_{j})= B_{j}(x),\qquad C(x,y_{j})= C_{j}(x)\quad \text{etc.}\end{gathered}
(8.62)

На прямых, описываемых уравнениями y_{j}=\text{const},~ j=1,\ldots,J-1, получаем


\begin{gathered} A_{j}(x)\widehat{u}\,''_{j}(x)+ B_{j}(x)\frac{\widehat{u}\,'_{j+1}(x)-\widehat{u}\,'_{j-1}(x)}{h}+ C_{j}(x) \frac{\widehat{u}_{j+1}(x)-2\,\widehat{u}_{j}(x)+ \widehat{u}_{j-1}(x)}{h^2}+ D_{j}(x) \widehat{u}\,'_{j}(x)+\\[2pt] E_{j}(x) \frac{\widehat{u}_{j+1}(x)-\widehat{u}_{j-1}(x)}{2h}+ F_{j}(x) \widehat{u}_{j}(x)= G_{j}(x),\quad j=1,\ldots,J-1;~ 0<x<L. \end{gathered}

Здесь \widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x) — приближенное решение задачи на прямыхy=y_1=\text{const},\ldots, y=y_{J-1}= \text{const}, A_{j}(x),\ldots, G_{j}(x) — коэффициенты уравнения (8.59) на этих прямых. В силу краевых условий находим


\begin{gathered}\widehat{u}_{0}(x)= u(x,0)= \varphi(x),\qquad \widehat{u}_{J}(x)= u(x,M)= \psi(x),\\[4pt] \widehat{u}_{j}(0)= \Phi(y_{j})= \Phi(j\,h),\quad \widehat{u}_{j}(L)= \Psi(y_{j})= \Psi(j\,h),\quad j=1,\ldots,J-1. \end{gathered}
(8.63)

Таким образом, от линейного дифференциального уравнения в частных производных (8.59) и краевыми условиями (8.60) переходим к дифференциально-разностной краевой задаче (8.62),(8.63) для системы из (n-1) обыкновенных дифференциальных уравнений (8.62) относительно (n-1) неизвестных функций \widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x). В результате ее решения можно найти законы изменения искомой функции u(x,y) вдоль прямых y_{j}=\text{const}, j=0,1,\ldots,J, то есть \widehat{u}(x,y_{j}),~ j=0,1,\ldots,J. Значения функции в промежуточных точках по у может быть найдено с помощью интерполяции. Система (8.62),(8.63) называется системой уравнений метода прямых.


Замечания

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

2. Метод прямых удобнее применять в том случае, когда коэффициенты в уравнении (8.59) не зависят от x. Тогда система (8.62) будет системой с постоянными коэффициентами.

3. Если расчетная область отлична от прямоугольника, при реализации метода возникают трудности, связанные с тем, что функции u_{j}(x) определены на различных множествах.

4. При решении нестационарных уравнений дискретизация проводится по пространственным переменным, а время t остается непрерывной переменной.

5. На практике система дифференциальных уравнений (8.62) из-за большого числа J должна решаться численно, поэтому потребуется вводить шаг не только по y.

6. Если изменение искомой функции внутри множества \Omega в направлении y не велико, то можно ограничиться одной или двумя прямыми.


Пример 8.9. Составить уравнения метода прямых для решения одномерной задачи теплопроводности (где \psi(x) — заданная функция)


\begin{aligned}& u_{t}=a^2\cdot u_{xx},\quad 0<x<1,\quad 0<t<1,\\ & u(x,0)= \psi(x),\quad 0 \leqslant x \leqslant 1,\\ & u(0,t)= u(1,t)=0,\quad 0<t \leqslant 1. \end{aligned}


Решение, формула (8.64)
Прямые, параллельные оси

Для решения задачи выберем на отрезке [0,1] точки x_{i}= i\,h,~ i=0,1,\ldots, I;~ I\,h=1 и проведем через них прямые, параллельные оси Ot (рис. 8.16). На каждой прямой приближенно заменим уравнение теплопроводности обыкновенным дифференциальным уравнением для функции u(x_{i},t). Для этого заменим частную производную u_{xx} по формуле


\left.{\frac{\partial^2 \widehat{u}}{\partial x^2}}\right|_{x=x_{i}}= \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}

и введем обозначения \widehat{u}(x_{i},t)= \widehat{u}_{i}(t),~ \frac{\partial \widehat{u}(x_{i}, t)}{\partial t}= \widehat{u}\,'_{i}(t).


Тогда на прямых x_{i}=\text{const},~ i=1,\ldots,I-1, получаем систему и


\begin{gathered}\widehat{u}\,'_{i}(t)= a^2\cdot \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}\,,\quad 1 \leqslant i \leqslant I-1,\\[4pt] \begin{aligned}& \widehat{u}_{0}(t)= u(0,t)=0, & & 0<t \leqslant 1,\\ & \widehat{u}_{I}(t)= u(1,t)=0, & & 0<t \leqslant 1 \end{aligned} \end{gathered}
(8.64)

с начальными условиями \widehat{u}_{i}(0)= \psi(x_{i})= \psi(i\,h),~ i=1,2,\ldots,I-1. Полученную систему дифференциальных уравнений можно записать в векторно-матричной форме:


\frac{d}{dt}\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1} \end{pmatrix}= \frac{a^2}{h^2}\! \begin{pmatrix}-1&1&\cdots& 0&0\\ 1&-2&1&\cdots&0\\ \vdots& \vdots& \vdots& \vdots& \vdots\\ \vdots& \vdots& 1&-2&1\\ 0&0&\cdots& 1&-2 \end{pmatrix}\!\cdot\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1}\end{pmatrix}\!.

Задача имеет аналитическое решение


\widehat{u}(t)= H^{-1}\cdot\! \begin{pmatrix}\exp \dfrac{\lambda_1a^2t}{h^2}&0&0\\ 0&\ddots&0\\ 0&0& \exp \dfrac{\lambda_{I-1}a^2t}{h^2} \end{pmatrix}\!\cdot H\,u(0),

где u(0)= \begin{pmatrix}\psi(h)\\ \psi(2h)\\ \vdots\\ \psi \bigl((I-1)h\bigr)\end{pmatrix}\!,~ \lambda_{i}=-\left(2\cos \frac{\pi\,i}{n}+ 2\right)\!,~ i=1,\ldots,I-1 — собственные значения матрицы системы, элементы матрицы H (собственные векторы матрицы A) равны H_{ik}= \sin \frac{\pi\,ik}{n}\,, i=1,\ldots,I-1, k=1,\ldots,I-1.


Заметим, что при численном решении системы обыкновенных дифференциальных уравнений при увеличении числа I (улучшении точности) коэффициент жесткости S системы увеличивается. При больших значениях n справедливо S= \frac{\max_{i}|\lambda_{i}|}{\min_{i}|\lambda_{i}|}\cong \frac{4I^2}{\pi^2}.




Метод характеристик


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


Рассмотрим основные идеи метода на примере задачи Коши для квазилинейной системы дифференциальных уравнений m-го порядка с двумя независимыми переменными x,\,t\colon


\begin{aligned}&\begin{aligned}&a_{11} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{1m} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{11} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{1m} \frac{\partial u_{m}(x,t)}{\partial x}= f_1(x,t,u),\\[-2pt] &\quad\vdots\\[-4pt] &a_{m1} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{mm} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{m1} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{mm} \frac{\partial u_{m}(x,t)}{\partial x}= f_m(x,t,u),\end{aligned}\\[4pt] &\begin{aligned}& u_1(x,0)= \varphi_1(x), & & x\in \mathbb{R},\\[-2pt] & \quad\vdots\\[-4pt] & u_m(x,0)= \varphi_m(x), & & x\in \mathbb{R}, \end{aligned} \end{aligned}

где a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ i=1,\ldots,m;~ j=1,\ldots,m;~ D=R\times (0,T),~ T>0 — заданный момент времени. Перепишем систему в векторно-матричном виде


\begin{gathered}A(x,t,u)\frac{\partial u(x,t)}{\partial t}+ B(x,t,u)\frac{\partial u_{1}(x,t)}{\partial x}= f(x,t,u),\quad (x,t)\in D,\\[2pt] u(x,0)= \varphi(x),\quad x\in \mathbb{R}, \end{gathered}
(8.65)

где A(x,t,u)= \bigl(a_{ij}(x,t,u)\bigr),~ B(x,t,u)= \bigl(b_{ij}(x,t,u)\bigr) — матрицы размеров (m\times m),~ f(x,t,u)= \begin{pmatrix}f_1(x,t,u)\\\vdots\\ f_m(x,t,u) \end{pmatrix} — вектор-функция, u(x,t)= \begin{pmatrix}u_1(x,t)\\ \vdots\\ u_m(x,t) \end{pmatrix} — искомая вектор-функция правых частей, \varphi(x)= \begin{pmatrix} \varphi_1(x)\\ \vdots\\ \varphi_m(x) \end{pmatrix} — заданная вектор-функция.


Пусть известно, что система (8.65) имеет гладкое решение в области D и на некоторой кривой \gamma оно задано и равно \Bigl.{u(x,t)}\Bigr|_{\gamma} (рис. 8.17). Рассмотрим задачу о нахождении решения в окрестности кривой \gamma.


График кривой, проходящей через точку

Выберем на кривой \gamma произвольную точку (x_0,t_0), а вектор бесконечно малого смещения вдоль этой кривой из точки (x_0,t_0) обозначим (dx,dt). Присоединим к уравнениям (8.65) выражение для полного дифференциала (для упрощения записи опустим аргументы):


A\,\frac{\partial u}{\partial t}+ B\,\frac{\partial u}{\partial x}=f,\qquad du= \frac{\partial u}{\partial x}\,dx+ \frac{\partial u}{\partial t}\,dt\,.

Умножим первое уравнение на dt, а из второго выразим \frac{\partial u}{\partial t}\,dt\colon


A\,\frac{\partial u}{\partial t}\,dt+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt,\qquad \frac{\partial u}{\partial t}\,dt= du-\frac{\partial u}{\partial x}\,dx\,.

Подставляя второе уравнение в первое, имеем A\,du-A\,\frac{\partial u}{\partial x}\,dx+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt или


(B\,dt-A\,dx)\cdot \frac{\partial u}{\partial x}= f\,dt-A\,du\,.
(8.66)

Если определитель матрицы (B\,dt-A\,dx) отличен от нуля в каждой точке кривой \gamma, то частная производная \frac{\partial u}{\partial x} и, следовательно, частная производная \frac{\partial u}{\partial t} определяются однозначно. В этом случае можно найти решение в окрестности кривой \gamma\colon\, \Bigl.{u(x,t)= u(x,t)}\Bigr|_{\gamma}+ du, что соответствует разложению искомого решения по формуле Тейлора до членов первого порядка. Если определитель матрицы (B\,dt-A\,dx) в каждой точке кривой \gamma равен нулю, то в силу предположения о существовании решения производные \frac{\partial u}{\partial x},\, \frac{\partial u}{\partial t} и, следовательно, функция u(x,t), находятся неоднозначно. В этом случае кривая у называется характеристикой. Приравнивая определитель матрицы (B\,dt-A\,dx) к нулю и предполагая, что dt\ne0 вдоль кривой \gamma, получаем соотношение


\left|B-A\cdot \frac{\partial u}{\partial t}\right|=0\quad \text{ili}\quad |B-A \lambda|=0,
(8.67)

где \lambda(x,t,u)= \frac{\partial u}{\partial t}. Линии, задаваемые дифференциалами смещения dx,\,dt, вдоль которых справедливо (8.67), являются характеристиками. Анализ разрешимости уравнения (8.67) позволяет ввести классификацию систем вида (8.65).


Если уравнение (8.67) имеет m различных действительных корней, то система (8.65) называется гиперболической; если их число меньше m, то система называется параболической; если уравнение не имеет действительных корней, то — эллиптической. В дальнейшем будем рассматривать только гиперболические системы, т.е. имеющие m различных характеристик:


\left(\frac{dx}{dt}\right)_{i}= \lambda_{i}(x,t,u),\quad i=1,2,\ldots,m,
(8.68)

где функции \lambda_{i}(x,t,u) задают характеристические направления в каждой точке кривой \gamma.


Вдоль характеристик решение системы должно удовлетворять характеристическим соотношениям (соотношениям совместности), которые получаются в результате приравнивания к нулю определителя m-го порядка расширенной матрицы, т.е.


\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr).
(8.69)

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


Этот способ получения характеристических соотношений основан на применении теоремы Кронекера-Капелли: решение системы (8.66) существует тогда и только тогда, когда ранг матрицы (B\,dt-A\,dx) этой системы равен рангу расширенной матрицы. Предполагается, что решение существует и определитель матрицы системы (m-го порядка) вдоль характеристик равен нулю. Поэтому должен быть равен нулю определитель расширенной матрицы, составленный из произвольных m столбцов.


В результате описан принцип получения характеристических уравнений (уравнений характеристик и характеристических соотношений). Их вывод составляет первый этап метода характеристик. Основные идеи первого, а также второго этапа, связанного с конструированием численной схемы, рассмотрим на примере линейной задачи Коши для одномерного волнового уравнения, описывающего распространение плоских звуковых волн в покоящейся среде.


Пример 8.10. Получить соотношения метода характеристик для решения линейной задачи Коши


\frac{\partial^2 \varphi(x,t)}{\partial t^2} = \frac{\partial^2 \varphi(x,t)}{\partial x^2}\,,\quad (x,t)\in \mathbb{R}\times (0;T),\quad \varphi(x,0)=x^2+x,\quad \varphi_t(x,0)=\cos x,\quad x\in \mathbb{R}.

Решение, формулы (8.70)-(8.77)

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


Для применения метода характеристик сначала сведем задачу к виду (8.65), т.е. к системе двух уравнений первого порядка. Для этого введем новые искомые функции


\nu(x,t)= \frac{\partial \varphi(x,t)}{\partial t}\,\qquad p(x,t)= \frac{\partial \varphi(x,t)}{\partial x}\,.

Здесь \nu(x,t),\,p(x,t) — отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн, \varphi(x,t) — потенциал функции \nu(x,t). Тогда уравнение можно представить в форме


\frac{\partial }{\partial t}\! \left(\frac{\partial \varphi(x,t)}{\partial t}\right)= \frac{\partial }{\partial x}\! \left(\frac{\partial \varphi(x,t)}{\partial x}\right)\quad \text{ili}\quad \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0.

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


\frac{\partial^2 \varphi(x,t)}{\partial x\,\partial t}= \frac{\partial^2 \varphi(x,t)}{\partial t\,\partial x}\quad \text{ili}\quad \frac{\partial p(x,t)}{\partial t}= \frac{\partial \nu(x,t)}{\partial x}\,.

В результате от исходной задачи для волнового уравнения переходим к задаче для системы двух дифференциальных уравнений с начальными условиями относительно новых искомых функций:


\left\{\begin{aligned}& \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0, & & \nu(x,0)=\cos x,\\ & \frac{\partial p(x,t)}{\partial t}-\frac{\partial \nu(x,t)}{\partial x}=0, & & p(x,0)=2x+1. \end{aligned}\right.

Запишем систему в виде (8.65)


A\cdot \frac{\partial u}{\partial t}+ B\cdot \frac{\partial u}{\partial x}=0,\quad u(x,0)= \begin{pmatrix}\cos x\\ 2x+1\end{pmatrix}\!,
(8.70)

где A=\begin{pmatrix}1&0\\0&1\end{pmatrix}\!,~ B=\begin{pmatrix}0&-1\\-1&0\end{pmatrix}\!,~ u=\begin{pmatrix}\nu\\p \end{pmatrix}\!,~ m=2. Для решения задачи (8.70) применим метод характеристик, содержащий два этапа.


Первый этап. Нахождение характеристик и вывод характеристических соотношений. Для этого решим уравнение (8.67), записанное для (8.70):


|B-A\cdot \lambda|= \begin{vmatrix}-\lambda&-1\\-1&-\lambda \end{vmatrix}=0.

Получаем \lambda^2-1=0,~ \lambda_{1,2}=\pm1, т.е. два направления характеристик:


{\left(\frac{dx}{dt}\right)\!}_{1}= \lambda_1=1;\quad {\left(\frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1.
(8.71)

Интегрируя, найдем общие решения полученных дифференциальных уравнений. Для первого уравнения x=t+\text{const}, а для второго x=-t+\text{const}. Им соответствуют два семейства характеристик (в данной задаче прямых линий). Их принято обозначать буквой C^{+} с верхним индексом "+" для семейства характеристик, имеющих меньший угол с положительным направлением оси Ох по сравнению с соответствующим углом для характеристик семейства C^{-}. В данной задаче характеристики семейства C^{+} имеют положительный наклон, а характеристики семейства C^{-} — отрицательный наклон (рис. 8.18). Условно их также называют характеристиками первого и второго семейств соответственно.


Семейства наклонных прямых

Для получения характеристических соотношений запишем расширенную матрицу (в решаемой задаче f=\begin{pmatrix}0\\0 \end{pmatrix}):


\bigl(B\,dt-A\,dx,\, f\,dt-A,\,du\bigr)= \left(\underbrace{\begin{pmatrix}0&-dt\\-dt&0 \end{pmatrix}}_{B\,dt}-\underbrace{\begin{pmatrix}dx&0\\ 0&dx \end{pmatrix}}_{A\,dx},-\underbrace{\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!\! \begin{pmatrix}d\nu\\ dp\end{pmatrix}}_{A\,du} \right)= \begin{pmatrix}-dx&-dt-d\nu\\-dt&-dx&-dp \end{pmatrix}\!.

Приравнивая определитель матрицы, составленной из последних двух столбцов, нулю, получаем


\begin{vmatrix}-dt&-d\nu\\-dx&-dp \end{vmatrix}= dp\cdot dt-d\nu\cdot dx=0.

Поделив на dt (предполагается, что dt\ne 0), имеем dp-d\nu\, \frac{dx}{dt}=0. Так как на прямых семейства C^{+} выполняется равенство {\left( \frac{dx}{dt} \right)\!}_{1}= \lambda_1=1, то соответствующие характеристическое соотношение имеет вид


dp-d\nu=0.
(8.72)

Аналогично на прямых семейства C^{-} справедливо {\left( \frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1 и


dp+ d\nu=0.
(8.73)

Обыкновенные дифференциальные уравнения (8.71)–(8.73) образуют систему уравнений метода характеристик. Полученные характеристические соотношения (8.72), (8.73) в данной задаче могут быть проинтегрированы. В результате получаются конечные соотношения, называемые инвариантами Римана. Далее они будут использованы при построении численной схемы


p-\nu=\text{const} (вдоль семейства C^{+}),

p+\nu=\text{const} (вдоль семейства C^{-}),

Результатом первого этапа является получение уравнений характеристик и характеристических соотношений (следствий из них):


\begin{aligned}& {\left( \frac{dx}{dt} \right)\!}_{1}=1 & & (x=t+\text{const}); & & p-\nu=\text{const},\\ & {\left( \frac{dx}{dt} \right)\!}_{2}=-1 & & (x=-t+\text{const}); & & p+\nu=\text{const}.\end{aligned}
(8.74)

К системе (8.74) следует добавить начальные условия \nu(x,0)=\cos x,~ p(x,0)=2x+1.


Второй этап. Конструирование численной схемы и вычисление значений искомых функций.


1. Для простоты изложения выделим на оси Ох отрезок [0;l] и на нем зададим три точки 0;0,5;1 (рис.8.19, а). Они соответствуют нулевому временному слою. Значения функций p(x,t),\, \nu(x,t) в этих точках будем обозначать p_{i}^0,\,\nu_{i}^0,~ i=0,1,2, а сами точки (x_{i},t^n) парами индексов (i,n), где i — номер точки по переменной x, а n — номер временного слоя. Для n=0 с учетом начальных условий \nu(x,0)=\cos x,~ p(x,0)=2x+1 получаем


\begin{aligned}& \nu_{0}^{0}= \nu(0;0)=1; & & \nu_{1}^{0}= \nu(0,\!5;0)=0,\!8776; & & \nu_{2}^{0}= \nu(1;0)=0,\!5403;\\[2pt] & p_{0}^{0}= p(0;0)=1; & & p_{1}^{0}= \nu(0,\!5;0)=2; & & p_{2}^{0}= p(1;0)=3. \end{aligned}

Конструкция численной схемы

2. Выведем разностные уравнения, определяющие положение точек на следующем временном слое и значения искомых функций в них. Система (8.74) заменяется системой разностных уравнений. Численное решение задачи строится на характеристической сетке, которая формируется одновременно с процессом нахождения решения. При этом обычно рассматривается трехточечный угловой шаблон с узлами "1","2","3", составленный отрезками характеристик (из семейств C^{+} и C^{-}), пересекающихся в рассчитываемой точке "3" (рис. 8.19,б). Назовем этот шаблон базовым. Предполагается, что положение узлов "1" , "2" и значения решения в узлах известны из начальных условий или рассчитаны на предыдущем шаге. Положение узла "3" и значение решения в нем подлежат определению.


Аппроксимируя производную по формулам


{\left(\frac{dx}{dt}\right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1} (вдоль характеристики из семейства C^{+} с \lambda_1=1),

{\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2} (вдоль характеристики из семейства C^{-} с \lambda_2=-1),

получаем

x_3-x_1=t_3-t_1,\qquad x_3-x_2=t_2-t_3.
(8.75)

Условие p-\nu=\text{const} вдоль характеристик из семейств C^{+} записывается для лежащих на ней узлов "1" и "3":


p_3-\nu_3=p_1-\nu_1,
(8.76)

а условие p-\nu=\text{const} вдоль характеристики из семейства C^{-} записывается для лежащих на ней узлов "2" и "3":


p_3+\nu_3= p_2+\nu_2.
(8.77)

Система (8.75)–(8.77) является замкнутой относительно искомых значений x_3,\,t_3,\,p_3,\,\nu_3 при известных x_1, t_1, x_2, t_2, p_1, \nu_1, p_2, \nu_2. Формулы (8.75)–(8.77) являются базовыми для расчета всех точек, принадлежащих области D. После совмещения узлов "1" и "2" базового шаблона с выбранными очередными узлами сетки в (8.75)–(8.77) подставляются конкретные величины в соответствии с принятыми обозначениями.


Чертеж углового шаблона

3. Производится расчет первого временного слоя. В решаемой задаче на нулевом временном слое угловой шаблон может занимать два положения (рис. 8.20, а,б). Для точки с номерами (0,1) первого рассчитываемого слоя узлы "1" и "2" углового шаблона совмещаются с узлами (0,0) и (1,0) (рис. 8,20,а). Поэтому система (8.75)–(8.77) имеет вид


\begin{aligned} & x_{0,1}-x_{0,0}= t_{0,1}-t_{0,0}, & & x_{0,1}-x_{1,0}= t_{1,0}-t_{0,1},\\[2pt] & p_0^1-\nu_0^1= p_0^0-\nu_0^0, & & p_0^1+ \nu_0^1= p_1^0+ \nu_1^0. \end{aligned}

Подставляя t_{0,0}=0,~ t_{1,0}=0,~ x_{0,0}=0,~ x_{1,0}=0,\!5 и известные значения p_0^0, p_1^0, \nu_0^0, \nu_1^0 на нулевом слое (см. п.1), получаем


\begin{gathered}x_{0,1}-0=t_{0,1}-0;\qquad x_{0,1}-0,\!5=0-t_{0,1};\\[2pt] p_0^1-\nu_0^1=1-1=0;\qquad p_0^1+\nu_0^1=2+ 0,\!8776= 2,\!8776. \end{gathered}

Решение этой системы: x_{0,1}=0,\!25;~ t_{0,1}=0,\!25;~ p_0^1=1,\!4388;~ \nu_0^1=1,\!4388.


Для второй точки первого слоя (рис. 8.20,б) с номерами (l;l) узлы "1" и "2" углового шаблона совмещаются с узлами (l;0) и (2;0). Тогда система (8.75)–(8.77) имеет вид:


\begin{gathered}x_{1,1}-x_{1,0}=t_{1,1}-t_{1,0};\qquad x_{1,1}-x_{2,0}= t_{2,0}-t_{1,1};\\[2pt] p_1^1-\nu_1^1= p_1^0-\nu_1^0;\qquad p_1^1+ \nu_1^1= p_2^0+ \nu_2^0. \end{gathered}

Подставляя t_{1,0}=0,~ t_{2,0}=0,~ x_{1,0}=0,\!5,~ x_{2,0}=1 и известные значения p_1^0,\nu_1^0,o_2^0,\nu_2^0 на нулевом слое (см. п.1), получаем


\begin{gathered}x_{1,1}-0,\!5=t_{11}-0;\qquad x_{1,1}-1=0-t_{1,1};\\[2pt] p_1^1-\nu_1^1=2-0,\!8776=1,\!1224;\qquad p_1^1+\nu_1^1=3+0,\!5403= 3,\!5403. \end{gathered}

Отсюда x_{1,1}= 0,\!75;~ t_{1,1}= 0,\!25; p_1^1= 2,\!33135;~ \nu_1^1= 1,\!20895.


4. Для расчета точки (0,2) следующего временного слоя узлы "1" и "2" углового шаблона совмещаются с узлами (0,1) и (1,1) (рис.8.20,в). Тогда система (8.75)–(8.77) принимает вид:


\begin{gathered}x_{0,2}-x_{0,1}=t_{0,2}-t_{0,1};\qquad x_{0,2}-x_{1,1}=t_{1,1}-t_{0,2};\\[2pt] p_0^2-\nu_0^2= p_0^1-\nu_0^1;\qquad p_0^2+\nu_0^2= p_1^1+ \nu_1^1. \end{gathered}

Подставляя значения с предыдущего слоя, имеем


\begin{gathered}x_{0,2}-0,\!25=t_{0,2}-0,\!25;\qquad x_{0,2}-0,\!75= 0,\!25-t_{0,2};\\[2pt] p_0^2-\nu_0^2= 1,\!4388-1,\!4388;\qquad p_0^2+ \nu_0^2=2,\!33135+ 1,\!20895. \end{gathered}

Отсюда x_{0,2}= 0,\!5;~ t_{0,2}= 0,\!5;~ p_0^2=\nu_0^2= 1,\!77015.


Очевидно (см. рис. 8.19,а), количество узлов на каждом последующем слое уменьшается на единицу. Расчет продолжается до замыкания области, так как на последнем слое останется один узел. В данном примере этим узлом является точка (0,2).


Замечание. Рассмотренная в примере 8.10 задача является всего лишь модельной, выбранной для пояснения основных принципов конструирования метода характеристик и его численной реализации. В вычислительной практике, как правило, ставятся нелинейные задачи. Проиллюстрируем некоторые особенности их решения на следующем примере.


Пример 8.11. Получить соотношения метода характеристик для решения нелинейной задачи Коши


\begin{gathered}\begin{aligned}&a_{11}\frac{\partial u_1(x,t)}{\partial t}+ a_{12}\frac{\partial u_2(x,t)}{\partial t}+ b_{11}\frac{\partial u_1(x,t)}{\partial x}+ b_{12}\frac{\partial u_2(x,t)}{\partial x}=f_1,\\ &a_{21}\frac{\partial u_1(x,t)}{\partial t}+ a_{22}\frac{\partial u_2(x,t)}{\partial t}+ b_{21}\frac{\partial u_1(x,t)}{\partial x}+ b_{22}\frac{\partial u_2(x,t)}{\partial x}=f_2,\end{aligned}\quad (x,t)\in D\\[4pt] u_1(x,0)= \varphi_1(x),\quad u_2(x,0)= \varphi_2(x),\quad x\in \mathbb{R},\end{gathered}

где a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ f_{i}=f_{i}(x,t,u),~ i=1,2;~ j=1,2;~ \varphi_{1}(x),~ \varphi_2(x) — заданные функции; D=\mathbb{R}\times (0,T),~ T>0 — заданный момент времени.


[reshenie=Решение, формулы (8.78)-(8.80)]

Поставленная задача является частным случаем задачи (8.65) при m=2. Применим метод характеристик, включающий два этапа.


Первый этап. Нахождение характеристик и характеристических соотношений. Уравнение (8.67) имеет вид


|B-A \lambda|= \begin{vmatrix}b_{11}-\lambda a_{11}& b_{12}-\lambda a_{12}\\ b_{21}-\lambda a_{21}& b_{22}-\lambda a_{22}\end{vmatrix}=0 или a\,\lambda^2+ b\,\lambda+c=0,

где a=a_{11}a_{12}-a_{12}a_{21},~ b=a_{21}b_{12}+ a_{12}b_{21}-a_{11}b_{22}-a_{22}b_{11},~ c=b_{11}b_{22}-b_{12}b_{21}. Решая полученное уравнение, находим выражения для двух характеристических направлений в каждой точке расчетной области:


\begin{aligned}&{\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}(x,t,u_1,u_2)= \frac{-b+ \sqrt{b^2-ac}}{2a}\quad (C^{+});\\[2pt] &{\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}(x,t,u_1,u_2)= \frac{-b-\sqrt{b^2-ac}}{2a}\quad (C^{-}). \end{aligned}
(8.78)

Сформируем расширенную матрицу


\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr)= \begin{pmatrix}b_{11}dt-a_{11}dx& b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{21}dt-a_{21}dx& b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{pmatrix}

Характеристические соотношения получаются в результате приравнивания нулю определителя второго порядка (m=2), составленного из второго и третьего столбцов расширенной матрицы:


\begin{vmatrix}b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{vmatrix}=0.

Раскрыв определитель и поделив на dt\ne0, получим


\left(e-a\,\frac{dx}{dt}\right)\!du_1+ g\,du_2+ h\,dt+ q\,dx=0,

где e=a_{11}b_{22}-a_{21}b_{12},~ g=a_{12}b_{22}-a_{22}b_{12},~ h=b_{12}f_2-b_{22}f_1,~ q=a_{22}f_1-a_{12}f_2. Последние два слагаемых перепишем в эквивалентной форме h\,dt+ q\,dx= \left(h+ q\,\frac{dx}{dt}\right)\!dt.


Поскольку для семейства характеристик C^{+} выполняется равенство {\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}, а для семейства C^{-} — равенство {\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}, то характеристические \соотношения имеют вид


\begin{aligned}& (e-a \lambda_1)\,du_1+ g\,du_2+ (h+q \lambda_1)\,dt=0\quad (C^{+});\\[2pt] & (e-a \lambda_2)\,du_1+ g\,du_2+ (h+q \lambda_2)\,dt=0\quad (C^{-}). \end{aligned}
(8.79)

В результате получена характеристическая система (8.78), (8.79) из четырех обыкновенных дифференциальных уравнений. Они позволяют найти значения t,x,u_1,u_2 вдоль характеристик.


Второй этап. Конструирование численной схемы. Как и в предыдущем примере, будем использовать трехточечный угловой шаблон (рис. 8.19, б). Предполагается, что в точках "1" и "2" известны значения t,x,u_1,u_2 из начальных данных или из расчета предыдущего слоя, т.е. t_1,x_1,u_{11},u_{21} для точки "1" и t_2,x_2,u_{12},u_{22} для точки "2". Здесь при обозначении значений функций u_1,u_2, а далее и при обозначении других функций в точках шаблона первый нижний индекс соответствует номеру самой функции, а второй индекс — номеру точки в шаблоне. Требуется найти значения t_3,x_3, u_{13}, u_{23} определяющие положение точки "3" и значения искомых функций в ней. Предполагается, что точки "1" и "2" близки друг к другу, поэтому заменим дифференциалы конечными приращениями, а производную \frac{dx}{dt} — конечноразностным отношением:


\begin{gathered}dt\cong \Delta t=t_3-t_1,\quad dx\cong \Delta x=x_3-x_1,\quad du_1\cong \Delta u_1=u_{13}-u_{11},\quad du_2\cong \Delta u_2=u_{23}-u_{21},\quad {\left(\frac{dx}{dt} \right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1}\quad (C^{+});\\[4pt] dt\cong \Delta t=t_3-t_2,\quad dx\cong \Delta x=x_3-x_2,\quad du_1\cong \Delta u_1=u_{13}-u_{12},\quad du_2\cong \Delta u_2=u_{23}-u_{22},\quad {\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2}\quad (C^{+}). \end{gathered}

В итоге уравнения характеристик являются уравнениями прямых, соединяющих точки "1" и "2" с точкой "3", а характеристическая система (8.78), (8.79) преобразуется к виду


\begin{gathered}\begin{aligned}& x_3-x_1= (\lambda_1)_1(t_3-t_1),\\ & x_3-x_2= (\lambda_2)_2(t_3-t_2),\end{aligned}\\[4pt] \begin{aligned}& (e-a \lambda_1)_1(u_{13}-u_{11})+ (g)_1(u_{23}-u_{21})+ (h+ q \lambda_1)_1 (t_3-t_1)=0,\\ & (e-a \lambda_2)_2(u_{13}-u_{12})+ (g)_2(u_{23}-u_{22})+ (h+ q \lambda_2)_2 (t_3-t_2)=0. \end{aligned}\end{gathered}
(8.80)

Здесь нижние индексы 1 и 2 указывают, что значения всех функций, входящих в соответствующие выражения, подсчитываются в точке "1" или "2". Решая систему (8.80), можно найти значения t_3,x_3,u_{13},u_{23}" определяющие положение точки "3" и значения искомых функций в ней. Сначала обычно решаются первые два уравнения и определяются значения t_3,x_3, а затем — последние два уравнения, откуда находятся u_{13},u_{23}.


Опишем теперь общую схему численного решения поставленной задачи Коши.


1. При t=0 на оси абсцисс задаются точки с некоторым шагом. В этих точках (на рис. 8.21 выделен фрагмент, содержащий точки 1-6) значения функций u_{1},u_{2} определяются начальными данными, т.е. функциями \varphi_1(x), \varphi_2(x).


Схема численного решения уравнения в частных производных

2. В соответствии с изложенной методикой последовательно рассчитываются координаты точек 7-11 и приближенные значения решения в них. При этом точки "1" и "2" углового шаблона последовательно совмещаются с точками 1 и 2, затем с точками 2 и 3, 3 и 4, 4 и 5, 5 и 6. Узлами выстраиваемой характеристической сетки служат точки пересечения характеристик (они образуют следующий временной слой).


3. Расчет слоев продолжается до достижения последней точки (точки 21) треугольной области.


Отметим, что дифференциальным уравнениям (8.78) соответствуют криволинейные характеристики, а точка "3" должна быть точкой пересечения этих криволинейных характеристик. При переходе от дифференциальных уравнений (8.78),(8.79) к разностным дифференциалы были заменены конечными разностями, а производные конечно-разностными отношениями. Вследствие этого криволинейные характеристики были аппроксимированы отрезками прямых. Решение системы (8.80) объявляется первым приближением, т.е. t_3^{(1)}= t_3,~ x_3^{(1)}= x_3, u_{13}^{(1)}= u_{13},~ u_{23}^{(1)}= u_{23}, и может возникнуть необходимость в уточнении координат точки "3" и значений искомых функций в этой точке. Ограничимся изложением одного из возможных способов уточнения. После нахождения первого приближения следует найти средние арифметические значения:


\begin{aligned}& \overline{\lambda}_{1}= \frac{(\lambda_1)_1+ (\lambda_1)_3}{2}\,,\quad \overline{\lambda}_{2}= \frac{(\lambda_2)_2+ (\lambda_2)_3}{2}\,,\quad K_1= \frac{(e-a \lambda_1)_1+ (e-a \lambda_1)_3}{2}\,,\quad K_2= \frac{(e-a \lambda_2)_2+ (e-a \lambda_2)_3}{2}\,,\\[4pt] & \overline{L}_{1}= \frac{(g)_1+ (g)_3}{2}\,,\quad \overline{L}_{2}= \frac{(g)_2+ (g)_3}{2}\,,\quad M_1= \frac{(h+ q \lambda_1)_1+ (h+ q \lambda_1)_3}{2}\,,\quad M_2= \frac{(h+ q \lambda_2)_2+ (h+ q \lambda_2)_3}{2}\,,\end{aligned}

где нижний индекс 3 указывает, что соответствующее выражение подсчитывается в найденной точке "3". Эти соотношения получаются, если при интегрировании дифференциальных уравнений применять правило трапеций. Затем решается система


\begin{gathered}x_{3}^{(2)}-x_{1}= \overline{\lambda}_{1}(t_{3}^{(2)}-t_{1}),\\ x_{3}^{(2)}-x_{2}= \overline{\lambda}_{2}(t_{3}^{(2)}-t_{2}),\\[4pt] K_1\,\bigl(u_{13}^{(2)}-u_{11}\bigr)+ L_1\,\bigl(u_{23}^{(2)}-u_{21}\bigr)+ M_1\,\bigl(t_{3}^{(2)}-t_{1}\bigr)=0,\\ K_2\,\bigl(u_{13}^{(2)}-u_{12}\bigr)+ L_2\,\bigl(u_{23}^{(2)}-u_{22}\bigr)+ M_2\,\bigl(t_{3}^{(2)}-t_{2}\bigr)=0 \end{gathered}

относительно нового (второго) приближения t_2^{(2)},\,x_3^{(2)},\, u_{13}^{(2)},\, u_{23}^{(2)}. Описанный процесс продолжают до тех пор, пока решения, полученные при двух последовательных приближениях, не будут совпадать с заданной точностью. Как правило, если расстояние между точками "1" и "2" невелико, достаточно сделать два уточнения.


Замечания

1. Приведенные примеры иллюстрируют важное свойство метода характеристик, отличающее его от ранее рассмотренных разностных методов. Оно заключается в том, что расчеты осуществляются вдоль характеристик — линий, по которым распространяется информация об исходных данных, что придает методу четкий физический смысл.

2. Описанный метод расчета, когда точки следующего слоя заранее неизвестны, а вычисления ведутся вдоль характеристик семейств C^{+} и C^{-}, выделяемых в процессе расчета, называется классическим методом характеристик. Существуют модификации этого метода, когда точки следующего слоя задаются заранее, а для нахождения значений искомых функций в них характеристики выпускаются в направлении к известному слою. Такие методы называются сеточно-характеристическими.




Метод Годунова


В механике сплошных сред при решении задач, описывающихся квазилинейными дифференциальными уравнениями в частных производных, в расчетных областях могут возникать поверхности сильного разрыва, к которым относятся ударные волны, центрированные волны разрежения, поверхности контактного разрыва. На этих поверхностях терпят разрыв как искомые функции, характеризующие исследуемые процессы, так и их производные. При решении таких задач используются два класса численных схем, в одном из которых с использованием законов сохранения явно выделяются поверхности разрыва, а в другом эти поверхности явно не выделяются, а рассчитываются "насквозь" вместе с областями непрерывности. К последнему классу методов относится метод С. К. Годунова. Идею его конструирования рассмотрим на примере решения смешанной краевой задачи для уравнений акустики.


Пример 8.12. Получить соотношения метода Годунова для решения задачи


\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0

с начальными условиями u(x,0)=\psi_1(x),~ p(x,0)=\psi_2(x),~ 0 \leqslant x \leqslant 1,
и краевыми условиями p(0,t)=\varphi_1(t),~ p(1,t)=\varphi_2(t),~ 0<t \leqslant T.

Решение, формулы (8.81)-(8.90)

Система двух дифференциальных уравнений описывает распространение плоских звуковых волн (малых возмущений) в покоящейся среде. Здесь u(x,t),\, p(x,t) — малые отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн; \rho_0 — плотность невозмущенной среды; a_0 — скорость звука.


Рассмотрим основные этапы конструирования метода Годунова.


1. Введем разбиение отрезка [0;1] с шагом h=1\!\!\not{\phantom{|}}\,n, где n — число промежутков. В результате получим узлы x_{i},~ i=0,1,\ldots,n, определяемые соотношением


x_{i+1}=x_{i}+h,\quad i=0,1,\ldots,n-1;\quad x_0=0,\quad x_n=1.

2. Заменим функции \psi_1(x),\, \psi_2(x) на начальном слое t=0 кусочно-постоянными (пример описанной замены для некоторой произвольной функции \psi_1(x) изображен на рис. 8.22) путем их интегрального осреднения по формулам


u_{1-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_1(x)\,dx\,,\quad p_{1-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_2(x)\,dx\,,\quad i=1,2,\ldots,n.
(8.81)

Здесь и далее все заданные или ранее рассчитанные функции снабжаются нижними полуцелыми индексами. Из рис. 8.22 видно, что полученная кусочно-постоянная функция имеет в узлах разрывы первого рода.


Кусочно-постоянные функции

3. Выберем шаг \tau по времени t и образуем следующий расчетный слой (t_1=\tau) с теми же узлами x_{i} (рис. 8.23). Таким образом, образуется совокупность прямоугольных ячеек. Предполагается, что искомые функции на горизонтальных и вертикальных границах образованных ячеек также являются постоянными.


Графики постоянных функций

На рассчитываемом слое при t=\tau приближенное решение задачи также представляется в виде кусочно-постоянной функции, имеющей разрывы в узлах сетки, совпадающих с узлами на предыдущем слое. Это означает, что во времени на слоях t=\text{const} сохраняется структура решения, изображенная на рис. 8.22. Значения искомых функций на рассчитываемом слое (верхних границах ячеек) обозначаются верхними полуцелыми индексами: u^{1-1\!\not{\phantom{|}}\,\,2},\, p^{1-1\!\not{\phantom{|}}\,\,2}, i=1,2,\ldots,n. На вертикальных границах ячеек значения функций обозначаются прописными буквами с нижними целыми индексами: U_{i},\,P_{i},~ i=0,\ldots,n.


При i=0,~ i=n значения P_0 и P_n определяются краевыми условиями по формулам, аналогичным (8.81):


P_0= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_1(t)\,dt\,,\quad P_n= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_2(t)\,dt\,.
(8.82)

Значения U_1,\ldots,U_{n-1},\, P_1,\ldots,P_{n-1} находятся в результате решения вспомогательной задачи о распаде произвольного разрыва, описанного в п.4, а значения U_0,\,U_n определяются в п.5.


4. Выберем две соседние ячейки с узлами x_{i-1},\,x_{i} и x_{i},\,x_{i+1} (эти ячейки выделены на рис. 8.23 жирными линиями). С целью нахождения значений U_{i},\,P_{i} на внутренней границе ячеек при x=x_{i} рассмотрим решение вспомогательной задачи, называемой задачей о распаде произвольного разрыва.


Дана система уравнений

\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0
(8.83)

с начальными условиями вида

u(x,0)= \left\{\begin{aligned}&u_{i-1\!\not{\phantom{|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &u_{i+1\!\not{\phantom{|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}, \end{aligned}\right.\qquad p(x,0)= \left\{\begin{aligned}&p_{i-1\!\not{\phantom{|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &p_{i+1\!\not{\phantom{|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}. \end{aligned}\right.

Требуется найти решение u(x,t),\,p(x,t) при t>0 и x=x_{i}, то есть u(x_{i},t)=U_{i},~ p(x_{i},t)=P_{i}.


Заметим, что в данной задаче начальные данные имеют разрыв первого рода в точке x_{i}. При t>0 они трансформируются, т.е. происходит распад разрыва, что обусловливает название рассматриваемой задачи. Для ее решения применим метод характеристик. Перепишем систему (8.83) в матричной форме


A\cdot \frac{\partial u(x,t)}{\partial t}+ B\cdot \frac{\partial u(x,t)}{\partial x}=0, где A=\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!,~~ B=\begin{pmatrix}0& 1\!\!\not{\phantom{|}}\,\rho_0\\ \rho_0a_0^2&0 \end{pmatrix}\!,~~ u=\begin{pmatrix}u\\p \end{pmatrix}.

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


|B-A\cdot \lambda|= \begin{vmatrix}-\lambda& 1\!\!\not{\phantom{|}}\,\rho_0\\ \rho_0a_0^2&-\lambda\end{vmatrix}=0. Отсюда {\left(\frac{dx}{dt}\right)\!}^{+}= \lambda^{+}= a_0,~~ {\left(\frac{dx}{dt}\right)\!}^{-}= \lambda^{-}=-a_0.

Интегрируя последние уравнения, имеем


x-a_0t=\text{const}~(C^{+}),\qquad x+a_0t=\text{const}~(C^{-}).

Если при t=0 выполняется условие x=x_{i}, то из последних соотношений выделяются две прямые, исходящие из точки x_{i} на начальном слое:


x-a_0t=x_{i}~(C^{+}),\qquad x+a_0t=x_{i}~(C^{-}).

Для получения характеристических соотношений запишем расширенную матрицу (8.69) (учитывая, что f\equiv 0):


\left(B\,dt-A\,dx,-A\! \begin{pmatrix}du\\dp\end{pmatrix}\right)= \begin{pmatrix}-dx&\dfrac{1}{\rho_0}\,dt&-dt\\ \rho_0a_0^2dt&-dx&-dp \end{pmatrix}\!.

Приравнивая определитель матрицы, составленной из двух последних столбцов, нулю, получаем


-\dfrac{1}{\rho_0}\,dp\,dt-dx\,du=0. Полагая dt\ne0, разделим уравнение на dt\colon\,-\dfrac{1}{\rho_0}\,dp-\frac{dx}{dt}\,du=0.

Подставляя вместо \frac{dx}{dt} сначала \lambda^{+}=a_0, а затем \lambda^{-}=-a_0, получаем


du+ \dfrac{1}{\rho_0\,a_0}\,dp=0\quad (C^{+}),\qquad du-\dfrac{1}{\rho_0\,a_0}\,dp= 0\quad (C^{-}).

Так как коэффициенты обоих уравнений постоянны, то их можно переписать в эквивалентной форме


d\! \left(u+ \dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{+}),\qquad d\! \left(u-\dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{-}).

Интегрируя и обозначая постоянные интегрирования I^{+},\,I^{-}, имеем


u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}\quad (C^{+}),\qquad u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}\quad (C^{-}).
(8.84)

Полученные соотношения называются инвариантами Римана. Они показывают, что величина u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+} остается постоянной вдоль семейства C^{+}, а величина u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-} остается постоянной вдоль семейства C^{-}.


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


Семества прямых, исходящие из узлов сетки

Рассмотрим две характеристики x-a_0t=x_{i-1} и x-a_0t=x_{i} из семейства C^{+}, исходящие из узлов x_{i-1},\,x_{i}, а также две характеристики x+a_0t=x_{i} и x+a_0t=x_{i+1} из семейства C^{-}, исходящие из узлов x_{i},\,x_{i+1} (рис. 8.24). В результате пересечения указанных характеристик образуются области I, II, III , первые две из которых представлены треугольниками, а третья — четырехугольником. Вертикальная граница (x_{i}= \text{const}), разделяющая ячейки, на которой ищется решение задачи, располагается в области III. Покажем, что на этой границе искомое решение постоянно. Для этого выберем на ней произвольную точку (назовем ее точка 3) , а решение в ней обозначим U_{i},\, P_{i}. В эту точку приходит характеристика семейства C^{+}, исходящая из точки 1 промежутка (x_{i-1}, x_{i}). Точке 1 соответствуют заданные значения u_{i-1\!\not{\phantom{|}}\,\,2},\, p_{i-1\!\not{\phantom{|}}\,\,2}. Вдоль этой характеристики остается постоянной величина u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}, то есть


U_{i}+ \frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i-1\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{i-1\!\not{\phantom{|}}\,\,2}.

Кроме того, в точку 3 приходит характеристика семейства C^{-}, исходящая из точки 2 промежутка (x_{i},x_{i+1}). Точке 2 соответствуют заданные значения и u_{i+ 1\!\not{\phantom{|}}\,\,2},\, p_{i+ 1\!\not{\phantom{|}}\,\,2}. Вдоль этой характеристики остается постоянной величина u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}, то есть


U_{i}-\frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i+ 1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0\,a_0}\cdot p_{i+ 1\!\not{\phantom{|}}\,\,2}.

Разрешая систему из двух последних уравнений относительно U_{i} и P_{i}, получаем


\begin{gathered}U_{i}= \frac{u_{i-1\!\not{\phantom{|}}\,\,2}+ u_{i+ 1\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{i+ 1\!\not{\phantom{|}}\,\,2}-p_{i+ 1\!\not{\phantom{|}}\,\,2}}{2\rho_0\,a_0}\,\\ P_{i}= \frac{p_{i-1\!\not{\phantom{|}}\,\,2}+ p_{i+ 1\!\not{\phantom{|}}\,\,2}}{2}+ \rho_0\,a_0\cdot \frac{u_{i+ 1\!\not{\phantom{|}}\,\,2}-u_{i-1\!\not{\phantom{|}}\,\,2}}{2}\,,\end{gathered}\quad i=1,2,\ldots,n-1.
(8.85)

Поскольку точка 3 была выбрана произвольно на вертикальной границе ячейки, то полученные выражения определяют решение в любой точке этой границы (очевидно, решение постоянно вдоль границы), за исключением самой точки x_{i} и точки пересечения характеристик, исходящих из точек x_{i-1},\, x_{i+1}. Таким образом, формулы (8.85) определяют "большие" величины U_{i},\, P_{i} на всех внутренних вертикальных границах ячеек.


Результаты решения данной вспомогательной задачи определяют шаг \tau по времени. Действительно, влияние каждого из разрывов, относящихся к соседним узлам x_{i-1},\,x_{i} на слое t=0, ограничивается характеристиками, исходящими из этих узлов. Физически это означает, что разрывы, возникшие при t=0, распространяются со скоростью звука a_0. Характеристики, исходящие из соседних узлов навстречу друг другу, пересекутся в точке с координатами t=\frac{h}{2a_0},~ x=\frac{x_{i-1}+ x_{i}}{2}. От этой встречи разрывов возникнут новые волны, которые через следующий отрезок времени длительностью \frac{h}{2a_0} достигнут вертикальных границ (x_{i-1}=\text{const},~ x_{i}=\text{const}) соседних ячеек. Именно до этого момента, т.е. в течение времени \frac{h}{a_0}= \frac{h}{2a_0}+ \frac{h}{2a_0}, значения U_{i},\, P_{i} будут сохраняться неизменными и определяться по формулам (8.85). Таким образом, величина шага по времени должна выбираться из условия \tau< \frac{h}{a_0}.


При i=0 для нахождения величины U_0 при известном P_0 используется инвариант Римана u-\frac{p}{\rho_0 a_0}=I^{-} вдоль семейства C^{-} (рис. 8.25). Поскольку на промежутке (0,x_1) заданы значения u_{1\!\not{\phantom{|}}\,\,2},\,p_{1\!\not{\phantom{|}}\,\,2}, то U_0-\frac{P_0}{\rho_0a_0}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{p_{1\!\not{\phantom{|}}\,\,2}}{\rho_0 a_0}. Отсюда получаем


U_0= \frac{1}{\rho_0a_0}\cdot P_0+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0 a_0}\cdot p_{1\!\not{\phantom{|}}\,\,2}.
(8.86)

Семейства прямых, образующих прямоугольные области

При i=n для нахождения величины U_n при известном P_n используется инвариант Римана u+\frac{p}{\rho_0a_0}=I^{+} вдоль семейства C^{+} (рис. 8.25). Поскольку на последнем промежутке (x_{n-1},1) заданы значения u_{n-1\!\not{\phantom{|}}\,\,2},\, p_{n-1\!\not{\phantom{|}}\,\,2}, то U_n+\frac{P_n}{\rho_0a_0}= u_{n-1\!\not{\phantom{|}}\,\,2}+ \frac{P_{n-1\!\not{\phantom{|}}\,\,2}}{\rho_0a_0}. Отсюда получаем


U_n=-\frac{1}{\rho_0a_0}\cdot P_n+ u_{n-1\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0a_0}\cdot p_{n-1\!\not{\phantom{|}}\,\,2}.
(8.87)

На рис. 8.25 показаны две крайние ячейки с тремя характеристиками, одна из которых (верхняя) ограничивает область влияния начальных данных, заданных на промежутках (x_0,x_1) и (x_{n-1},x_{n}) соответственно.


6. Перейдем к задаче расчета искомых величин на следующем временном слое. Для этого рассмотрим одну ячейку (рис. 8.26), для которой известны значения приближенного решения u_{i-1\!\not{\phantom{|}}\,\,2},\, p_{i-1\!\not{\phantom{|}}\,\,2} на нижней границе и U_{i-1},P_{i-1},U_{i},P_{i}, i=1,2,\ldots,n, на левой и правой вертикальных границах (см. п.4). Требуется определить значения p^{i-1\!\not{\phantom{|}}\,\,2},\, u^{i-1\!\not{\phantom{|}}\,\,2} на верхней границе ячейки.


Прямоугольная область и ячейка сетки

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


\frac{\partial u(x,t)}{\partial t}+ \frac{\partial}{\partial x}\! \left(\frac{p(x,t)}{\rho_0} \right)=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \frac{\partial}{\partial x}\bigl(\rho_0a_0^2u(x,t) \bigr)=0.
(8.88)

Запишем двойной интеграл от левой и правой частей уравнений по поверхности ячейки. Для его нахождения применим формулу Грина перехода от двойного интеграла к криволинейному интегралу второго рода:


\iint\limits_{S}\! \left(\frac{\partial \Phi_1}{\partial t}+ \frac{\partial \Phi_2}{\partial x}\right)\!dx\,dt= \oint\limits_{\gamma}\bigl(\Phi_1\,dx-\Phi_2\,dt\bigr),

где S — поверхность ячейки, а \gamma — ее граница, обходимая против часовой стрелки. Таким образом, вместо системы (8.88) получается система интегральных уравнений


\oint\limits_{\gamma}u\,dx-\frac{p}{\rho_0}\,dt=0,\qquad \oint\limits_{\gamma}p\,dx-\rho_0a_0^2u\,dt=0.

Проводя интегрирование, имеем


\begin{gathered}u_{i-1\!\not{\phantom{|}}\,\,2}h-\frac{P_{i}}{\rho_0}\cdot \tau-u^{i-1\!\not{\phantom{|}}\,\,2}h+ \frac{P_{i-1}}{\rho_0}\cdot \tau=0,\\[4pt] p_{i-1\!\not{\phantom{|}}\,\,2}h-\rho_0a_0^2U_{i}\tau-p^{i-1\!\not{\phantom{|}}\,\,2}h+ \rho_0a_0^2U_{i-1}\tau=0. \end{gathered}

Из этих уравнений получаем


\begin{gathered}u^{i-1\!\not{\phantom{|}}\,\,2}= u_{i-1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}\cdot (P_{i}-P_{i-1}),\\[4pt] p^{i-1\!\not{\phantom{|}}\,\,2}= p_{i-1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\cdot \rho_0a_0^2(U_{i}-U_{i-1}). \end{gathered}

Таким образом, получены формулы для расчета первого временного слоя (t=\tau). После нахождения значений u^{\frac{1}{2}},p^{\frac{1}{2}}, u^{\frac{3}{2}}, p^{\frac{3}{2}},\ldots, u^{n-\frac{1}{2}}, p^{n-\frac{1}{2}} из условия \tau<\frac{h}{a_0} выбирается очередной шаг и формируется новый временной слой (t_2=t_1+ \tau). Найденные значения используются в качестве опорных (они переобозначаются путем присвоения им соответствующих нижних индексов) для нахождения приближенного решения на следующем слое по формулам (8.85)–(8.87), (8.89). Заметим, что для вычисления P_0,\,P_n вместо (8.82) используются формулы интегрального осреднения, аналогичные (8.81). Для расчета, например, второго временного слоя они имеют вид


P_0= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_1(t)\,dt\,,\qquad P_n= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_2(t)\,dt\,.
(8.90)

Изложенная схема метода Годунова имеет первый порядок аппроксимации по \tau и h.


Замечание. Из процедуры интегрирования следует, что искомые значения u^{n-1\!\not{\phantom{|}}\,\,2},\, p^{n-1\!\not{\phantom{|}}\,\,2} соответствуют среднеинтегральным величинам


u^{n-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} u(x,\tau)\,dx\,,\qquad p^{n-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} p(x,\tau)\,dx\,.

Приведем пример расчета по описанной схеме для частного случая рассматриваемой задачи (где \rho_0=0,\!123;~ a_0=0,\!25):


\begin{gathered}\begin{aligned}& \frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\\ & \frac{\partial p(x,t)}{\partial t}+\rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0,\end{aligned}\qquad 0<x<1,\quad 0<t<1,\\[4pt] \begin{aligned}& u(x,0)= \psi_1(x)=\sin x, & & p(x,0)=\psi_2(x)= \rho_0x\sin x, & & 0 \leqslant x \leqslant 1,\\ & p(0,t)= \varphi_1(t)=0,\!1t+ \rho_0a_0^2\sin t, & & p(1,t)= \varphi_2(t)= 0,\!1t+ \rho_0a_0^2\sin t, & & 0<t \leqslant1, \end{aligned} \end{gathered}

Решение

А. Выберем число промежутков n=2. Тогда находится шаг h=\frac{1}{2}= 0,\!5 и вводятся узлы x_0=0;~ x_1=0,\!5;~ x_2=1. Зададим шаг \tau= 0,\!5 по времени из условия \tau< \frac{h}{a_0}= \frac{0,\!5}{0,\!25}=2. Сформируем две ячейки "A" и "B" первого временного слоя t_1=\tau= 0,\!5 (рис. 8.27).


Ячейки сетки в прямоугольной области

Б. Заменим функции \psi_1(x)=\sin x,~ \psi_2(x)=\rho_0x\sin x кусочно-постоянными по формулам (8.81). Для вычисления интегралов используются квадратурные формулы трапеций:


u_{i-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h}\cdot \frac{u_{i-1}+ u_{i}}{2}\cdot h=\frac{u_{i-1}+ u_{i}}{2}\,,\qquad p_{i-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h}\cdot \frac{p_{i-1}+ p_{i}}{2}\cdot h=\frac{p_{i-1}+ p_{i}}{2}\,,\quad i=1;2,

где u_{i}=\sin x_{i},~ p_{i}=\rho_0x_{i}\sin x_{i}. Поэтому


\begin{aligned}& u_{1\!\not{\phantom{|}}\,\,2}\cong \frac{0+0,\!479}{2}\cong 0,\!239; & & u_{3\!\not{\phantom{|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660;\\[2pt] <br />& p_{1\!\not{\phantom{|}}\,\,2}\cong \frac{0+0,\!0294}{2}\cong 0,\!239; & & p_{3\!\not{\phantom{|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660<br />\end{aligned}

В. Найдем значения P_0,\,P_2 по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций:


P_0\cong \frac{0+0,\!0537}{2}\cong 0,\!0269;\qquad P_2\cong \frac{0+0,\!0463}{2}\cong 0,\!0231.

Рассчитаем значение U_0 на левой границе ячейки "А" по формуле (8.86): U_0=\frac{P_0}{\rho_0a_0}+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{p_{1\!\not{\phantom{|}}\,\,2} }{\rho_0a_0}\cong 0,\!635 и значение U_2 на правой границе ячейки по формуле (8.87): U_2=-\frac{P_2}{\rho_0a_0}+ u_{3\!\not{\phantom{|}}\,\,2}+ \frac{p_{3\!\not{\phantom{|}}\,\,2} }{\rho_0a_0}\cong 1,\!476.


Г. Вычислим значения U_1,\,P_1 на правой границе ячейки "А" (левой границе ячейки "В" ). По формулам (8.85) получаем


U_1=\frac{u_{1\!\not{\phantom{|}}\,\,2}+ u_{3\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{|}}\,\,2}-p_{1\!\not{\phantom{|}}\,\,2}}{2\rho_0a_0}\cong-0,\!392;\qquad \frac{p_{1\!\not{\phantom{|}}\,\,2}+ p_{3\!\not{\phantom{|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{|}}\,\,2}-u_{1\!\not{\phantom{|}}\,\,2}}{2}\cong 0,\!0341.

Д. Вычислим значения u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2} на верхних границах ячеек:


\begin{aligned}& u^{1\!\not{\phantom{|}}\,\,2}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!1805; & & p^{1\!\not{\phantom{|}}\,\,2}= p_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_1-U_0)\cong 0,\!0226;\\[4pt] & u^{3\!\not{\phantom{|}}\,\,2}= u_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!749; & & p^{3\!\not{\phantom{|}}\,\,2}= p_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_2-U_1)\cong 0,\!0521. \end{aligned}

Таким образом, первый слой полностью рассчитан.


Е. Для расчета второго слоя (t_2=t_1+\tau) с шагом \tau=0,\!5 полученные значения u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2} переобозначим, для чего верхние индексы заменяются нижними:


u_{1\!\not{\phantom{|}}\,\,2}=0,\!1805;\quad p_{1\!\not{\phantom{|}}\,\,2}= 0,\!0226;\quad u_{3\!\not{\phantom{|}}\,\,2}=0,\!749;\quad p_{3\!\not{\phantom{|}}\,\,2}= 0,\!0521

Сформируем две ячейки "С" и "D" второго временного слоя (см. рис. 8.27).


Ж. Найдем значения P_0,\,P_2 по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций:


P_0\cong \frac{0,0537+ 0,1065}{2}\cong 0,\!0801;\qquad P_2\cong \frac{0,0463+ 0,00935}{2}\cong 0,\!0278.

Рассчитаем значение U_0 на левой границе ячейки "С" по формуле (8.86):


U_0= \frac{1}{\rho_0\,a_0}\cdot P_0+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0\, a_0}\cdot p_{1\!\not{\phantom{|}}\,\,2}\cong 2,\!0504

и значение U_2 на правой границе ячейки "D" по формуле (8.87):


U_2=-\frac{1}{\rho_0\,a_0}\cdot P_2+ u_{3\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{3\!\not{\phantom{|}}\,\,2}\cong 0,\!8651.

3. Вычислим значения U_1,\,P_1 на правой границе ячейки "С" (левой границе ячейки "D"). По формулам (8.85) получаем


U_1= \frac{u_{1\!\not{\phantom{|}}\,\,2}+ u_{3\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{|}}\,\,2}-p_{1\!\not{\phantom{|}}\,\,2}}{2\,\rho_0\,a_0}\cong-0,\!0150;\qquad P_1=\frac{p_{1\!\not{\phantom{|}}\,\,2}+ p_{3\!\not{\phantom{|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{|}}\,\,2}-u_{1\!\not{\phantom{|}}\,\,2}}{2}\cong 0,\!02861.

Вычислим значения u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2} на верхних границах ячеек "С" и "D":


\begin{aligned}& u^{1\!\not{\phantom{|}}\,\,2}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!599; & & p^{1\!\not{\phantom{|}}\,\,2}= p_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_1-U_0)\cong 0,\!0385;\\[2pt] & u^{3\!\not{\phantom{|}}\,\,2}= u_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!7556; & & p^{3\!\not{\phantom{|}}\,\,2}= p_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_2-U_1)\cong 0,\!04533. \end{aligned}

Таким образом, второй слой полностью рассчитан и найдено приближенное решение поставленной задачи.


Замечание. Здесь рассмотрен пример методики конструирования метода С.К. Годунова и ее реализации для относительно простой линейной задачи. В нелинейном случае метод существенно усложняется. Для более глубокого изучения проблемы авторы рекомендуют, где данный метод применяется для различных классов нелинейных задач и задач с тремя и четырьмя независимыми переменными.

Математический форум (помощь с решением задач, обсуждение вопросов по математике).
Кнопка "Поделиться"
Если заметили ошибку, опечатку или есть предложения, напишите в комментариях.

Часовой пояс: UTC + 3 часа [ Летнее время ]


Яндекс.Метрика

Copyright © 2010-2023 MathHelpPlanet.com. All rights reserved