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

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

Теоретический раздел
Часовой пояс: UTC + 3 часа [ Летнее время ]
MathHelpPlanet.com RSS-лента Математического форума

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


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

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


Рассмотрим задачи с тремя независимыми переменными x,y,t, где x,y — пространственные переменные, а t — время, t=0 и t=T — моменты начала и окончания процесса. В качестве примера расчетного множества, в котором ищется решение задачи, используем цилиндр (рис. 9.1). В начальный момент времени t=0 множество задано односвязной областью \Omega с кусочно-гладкой границей \Gamma.


Множество, задано односвязной областью

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


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


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


\frac{\partial u(x,y,t}{\partial t}+ c_1\cdot \frac{\partial u(x,y,t}{\partial x}+ c_2\cdot \frac{\partial u(x,y,t}{\partial y}= f(x,y,t),\quad (x,y,t)\in \Omega\times (0,T),
(9.1)

где c_1,\,c_2 — скорости переноса по осям Ox,\,Oy,~f(x,y,t) — заданная функция. Оно является модельным для математического описания многих физических процессов (нейтронных потоков в реакторе; теплопроводности в газах, обусловленной диффузией атомов и электронов, и др.).


Начально-краевая задача:
u_t+u_x+u_y=f,\quad (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi(x,y),\quad (x,y)\in \Omega\cup \Gamma, (начальное условие)
u(x,y,t)= \varphi(x,y,t),\quad (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное начальное условие (при t=0) и функциональное краевое условие. Здесь и далее \psi(x,y),~ \varphi(x,y,t) — заданные функции. Требуется найти функцию u(x,y,t), удовлетворяющую дифференциальному уравнению (9.1), начальному и краевому условиям.




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


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


\frac{\partial u(x,y,t)}{\partial t}= a^2\! \left(\frac{\partial^2 u(x,y,t)}{\partial x^2}+ \frac{\partial^2 u(x,y,t)}{\partial y^2}\right)\!,\quad (x,y,t)\in \Omega\times (0,T),
(9.2)

где {u} —температура, a^2 — коэффициент температуропроводности, a^2=\frac{k}{\rho\,c},~ k — коэффициент теплопроводности, c — удельная теплоемкость, \rho — плотность. Оно описывает, например, распределение температуры на плоской пластине при отсутствии тепловых источников.


Первая начально-краевая задача:
u_t=a^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
u(x,y,t)= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное начальное условие (при t=0) и функциональное краевое условие на границе \Gamma. Здесь и далее \psi(x,y),\, \varphi(x,y,t) — заданные функции.


Вторая начально-краевая задача:
u_t=a^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\frac{\partial u(x,y,t)}{\partial n}= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное начальное условие (при t=0) и дифференциальное краевое условие на границе \Gamma, задающее производную в направлении внешней нормали.


Третья начально-краевая задача:
u_t=a^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\alpha\,\frac{\partial u(x,y,t)}{\partial n}+ \beta\,u= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное начальное условие (при t=0) и функционально-дифференциальное краевое условие, задающее линейную комбинацию производной в направлении внешней нормали и функции на фанице \Gamma. Здесь \alpha,\, \beta — заданные числа (\alpha^2+ \beta^2>0) либо функции от t в более общем случае.


Во всех перечисленных задачах требуется найти функцию u(x,y,t), которая удовлетворяет дифференциальному уравнению (9.2) в области D= \Omega\times (0,T) и соответствующим условиям на ее границе.


Примером физической задачи, приводящей к первой начально-краевой задаче, может быть процесс изменения распределения температуры в плоской пластине, на границе которой поддерживается заданная температура, определенная функцией \varphi(x,y,t). В начальный момент времени при t=0 функцией \psi(x,y) задано начальное распределение температуры на поверхности пластины. Искомая функция u(x,y,t) характеризует температуру пластины в точке с координатами (x,y) в некоторый момент времени t.




Постановки задач для двумерного волнового уравнения


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


\frac{\partial^2 u(x,y,t)}{\partial t^2}= c^2\! \left(\frac{\partial^2 u(x,y,t)}{\partial x^2}+ \frac{\partial^2 u(x,y,t)}{\partial y^2}\right)\!,\quad (x,y,t)\in \Omega\times (0,T),
(9.3)

где u(x,y,t) — величина смещения точки (x,y) в момент времени t,~ c^2=\widetilde{T}\!\!\not{\phantom{|}}\,\rho,~ \widetilde{T} — проекция вектора натяжения на плоскость Oxy,~\rho — поверхностная плотность мембраны. Волновое уравнение описывает малые поперечные колебания мембраны.


Первая начально-краевая задача:
u_{tt}= c^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi_1(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\frac{\partial u(x,y,0)}{\partial t}= \psi_2(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
u(x,y,t)= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное и дифференциальное начальные условия (при t=0) и функциональное краевое условие. Здесь и далее \psi_1(x,y),\,\psi_2(x,y),\,\varphi(x,y,t) — заданные функции.


Вторая начально-краевая задача:
u_{tt}= c^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi_1(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\frac{\partial u(x,y,0)}{\partial t}= \psi_2(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\frac{\partial u(x,y,t)}{\partial n}= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное и дифференциальное начальные условия (при t=0) и дифференциальное краевое условие на границе \Gamma, задающее производную в направлении внешней нормали.


Третья начально-краевая задача:
u_{tt}= c^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),
u(x,y,0)= \psi_1(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\frac{\partial u(x,y,0)}{\partial t}= \psi_2(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)
\alpha\,\frac{\partial u(x,y,t)}{\partial n}+ \beta\,u= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)

содержит функциональное и дифференциальное начальные условия (при t=0) и функционально-дифференциальное краевое условие, задающее линейную комбинацию производной в направлении внешней нормали и функции на границе \Gamma. Здесь \alpha,\,\beta — заданные числа (\alpha^2+\beta^2>0) либо функции от t в более общем случае.


Во всех перечисленных задачах требуется найти функцию u(x,y,t), которая удовлетворяет дифференциальному уравнению (9.3) в области D= \Omega\times (0,T) и соответствующим условиям на ее границе.




Постановки задач для трехмерного уравнения Лапласа


Рассмотрим трехмерное уравнение Лапласа


\frac{\partial^2 u(x,y,t)}{\partial x^2}+ \frac{\partial^2 u(x,y,t)}{\partial y^2}+ \frac{\partial^2 u(x,y,t)}{\partial z^2}=0,\quad (x,y,t)\in \Omega,
(9.4)

или более общее уравнение Пуассона (где \Omega — некоторая односвязная область пространства Oxyz с кусочно-гладкой границей \Gamma,\, f(x,y,z) — заданная функция)


\frac{\partial^2 u(x,y,t)}{\partial x^2}+ \frac{\partial^2 u(x,y,t)}{\partial y^2}+ \frac{\partial^2 u(x,y,t)}{\partial z^2}= f(x,y,t),\quad (x,y,t)\in \Omega,

Поскольку в уравнении (9.4) отсутствуют производные по времени, задача является стационарной, следовательно, начальные условия в ней не задаются. На границе \Gamma области \Omega могут задаваться краевые (граничные) условия трех видов.


Первая краевая задача (задана Дирихле):
u_{xx}+u_{yy}+u_{zz}=0,~~ (x,y,z)\in \Omega,
\Bigl.{u}\Bigr|_{\Gamma}= \varphi(x,y,z),~~ (x,y,z)\in \Gamma, (краевое условие)

содержит функциональное краевое условие на границе области, \varphi(x,y,z) — заданная функция.


Вторая краевая задача (задача Неймана):
u_{xx}+u_{yy}+u_{zz}=0,~~ (x,y,z)\in \Omega,
\left.{\frac{\partial u}{\partial n}}\right|_{\Gamma}= \varphi(x,y,z),~~ (x,y,z)\in \Gamma, (краевое условие)

содержит дифференциальное краевое условие, в котором задается производная в направлении внешней нормали на границе \Gamma области \Omega.


Третья краевая задача:
u_{xx}+u_{yy}+u_{zz}=0,~~ (x,y,z)\in \Omega,
\left.{\alpha\,\frac{\partial u}{\partial n}}\right|_{\Gamma}+ \Bigl.{\beta\,u}\Bigr|_{\Gamma}= \varphi(x,y,z),~~ (x,y,z)\in \Gamma, (краевое условие)

содержит функционально-дифференциальное краевое условие, в котором задается линейная комбинация производной в направлении внешней нормали и функции u(x,y) на границе \Gamma области \Omega.


Во всех перечисленных задачах требуется найти функцию u(x,y), которая удовлетворяет дифференциальному уравнению (9.5) внутри области \Omega и соответствующим условиям на ее границе. В качестве примера задачи Дирихле можно привести задачу о нахождении стационарного распределения температуры внутри объемной области, если задана температура на ее границе.




Разностный метод


Принципы построения разностных схем изложены ранее. Для простоты изложения будем предполагать, что в двумерных начально-краевых задачах для уравнения переноса, уравнения теплопроводности и волнового уравнения множество D представляет собой параллелепипед длины L , ширины M и высоты T, ограниченный отрезками прямых, параллельных осям Ox,\,Oy и Ot в пространстве трех независимых переменных x,\,y и t (рис. 9.2). Поэтому зададим сетку D_h= \bigl\{x_{i},y_{j},t^n\bigr\}, где x_{i}=i\,h_x,~ 0 \leqslant i \leqslant I; y_{j}=j\,h_y,~ 0 \leqslant j \leqslant J; t^n=n\,\tau,~ 0 \leqslant n \leqslant N; N,\,J,\,I — целые положительные числа; h_x,\,h_y,\,\tau — величины шагов по пространственным координатам и времени (для простоты принимаются постоянными); I\,h_x=L, J\,h_y=M,~ N\,\tau=T. Такая сетка называется равномерной (регулярной).


Множество, задающие параллелепипед

Введем сеточную функцию u_h= \bigl\{u_{i,j}^{n}= u(x_{i},y_{j},t^n),~ 0 \leqslant i \leqslant I,~ 0 \leqslant j \leqslant J,~ 0 \leqslant n \leqslant N \bigr\}, которая является сеточным представлением решения исходной (дифференциальной) задачи или точным решением дифференциальной задачи в узлах сетки. Как правило, вычислить u_h не удается, поэтому находят другую сеточную функцию \widehat{u}_h\cong u_h, приближенно совпадающую с точным решением в узлах сетки. Она ищется как решение разностной схемы. Под разностной схемой, как и ранее, понимается совокупность разностных уравнений, аппроксимирующих основное дифференциальное уравнение во всех внутренних узлах сетки и дополнительные условия (начальные и краевые) — в граничных узлах.


Для аппроксимации частных производных можно использовать естественные обобщения формул (8.20)-(8.22), (8.26)-(8.32):


\begin{aligned}&\widehat{u}_{t}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\tau}\,, & &\widehat{u}_{t}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n}-u_{i,j}^{n-1}}{\tau}\,,\\[2pt] &\widehat{u}_{t}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n+1}-u_{i,j}^{n-1}}{2\tau}\,, & &\widehat{u}_{tt}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n+1}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{\tau^2}\, \end{aligned}
(9.5)

\begin{aligned}&\widehat{u}_{x}(x_{i},y_{j},t^n)= \frac{u_{i+1,j}^{n}-u_{i,j}^{n}}{h_x}\,, & &\widehat{u}_{x}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{h_x}\,,\\[2pt] &\widehat{u}_{x}(x_{i},y_{j},t^n)= \frac{u_{i+1,j}^{n}-u_{i-1,j}^{n}}{2h_x}\,, & &\widehat{u}_{xx}(x_{i},y_{j},t^n)= \frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{h_x^2}\,, \end{aligned}
(9.6)

\begin{aligned}&\widehat{u}_{y}(x_{i},y_{j},t^n)= \frac{u_{i,j+1}^{n}-u_{i,j}^{n}}{h_y}\,, & &\widehat{u}_{y}(x_{i},y_{j},t^n)= \frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{h_y}\,,\\[2pt] &\widehat{u}_{y}(x_{i},y_{j},t^n)= \frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2h_y}\,, & &\widehat{u}_{yy}(x_{i},y_{j},t^n)= \frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{h_y^2}\,. \end{aligned}
(9.7)

Заметим, что приведенные в (9.5)–(9.7) операторы, в знаменателе которых находятся \tau,\,h_x,\,h_y, имеют первый порядок аппроксимации по соответствующим шагам, а операторы с \tau^2,h_x^2,h_y^2,2\tau,2h_x,2h_y в знаменателе — второй порядок.


Пример 9.1. Построить разностную схему для задачи (где f(x,y,t),~ \psi(x,y),~ \varphi(y,t),~ \Phi(x,t) — заданные функции):
u_t+u_x+u_y= f(x,y,t),~~ (x,y,t)\in (0,L)\times (0,M)\times (0,T),
u(x,y,0)= \psi(x,y),~~ (x,y)\in [0,L]\times [0,M], (начальное условие)
u(0,y,t)=\varphi(y,t),~~ y\in [0,M],~ 0<t \leqslant T, (краевое условие)
u(x,0,t)=\Phi(x,t),~~ x\in [0,L],~ 0<t \leqslant T (краевое условие).

▼ Решение, формулы (9.8)-(9.9)

Решается первая начально-краевая задача для двумерного уравнения перс-носа (9.1), где c_1=c_2=t. Построим аналог абсолютно устойчивой схемы (8.40). Для этого используем сетку D_h и четырехточехный шаблон (рис. 9.3). Для аппроксимации производной u_t используем первую формулу в (9.5) и вторые формулы в (9.6) и (9.7), записанные на (n+1)-м временном слое (где f_{i,j}^{n}= f(x_{i},y_{j},t^{n}) — сеточное представление правой части):


\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n}}{\tau}+ \frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i-1,j}^{n+1}}{h_x}+ \frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j-1}^{n+1}}{h_y}= f_{i,j}^{n},\quad \begin{matrix}1 \leqslant i \leqslant I-1,~ 1 \leqslant j \leqslant J-1,\\ 1 \leqslant n \leqslant N-1.\end{matrix}
(9.8)

Четырехточехный шаблон

Начальные и краевые условия представляются в виде (где \varphi_{j}^n= \varphi(y_{j},t^{n}),~ \Phi_{i}^{n}= \Phi(x_{i},t^{n}),~ \psi_{i,j}=\psi(x_{i},y_{j}))


\begin{aligned}& \widehat{u}_{i,j}^{0}= \psi_{i,j}, & & 0 \leqslant i \leqslant I, & & 0 \leqslant j \leqslant J,\\[2pt] & \widehat{u}_{0,j}^{n}= \varphi_{j}^{n}, & & 0 \leqslant j \leqslant J, & & 1 \leqslant n \leqslant N,\\[2pt] & \widehat{u}_{i,0}^{n}= \Phi_{i}^{n}, & & 0 \leqslant i \leqslant I, & & 1 \leqslant n \leqslant N. \end{aligned}
(9.9)

Совокупность соотношений (9.8),(9.9) составляет неявную двухслойную разностную схему. Можно показать, что она абсолютно устойчивая и имеет первый порядок аппроксимации. Следовательно, она сходится с первым порядком точности. Последовательность вычислений по разностной схеме может быть различной. Достаточно, чтобы она соответствовала какому-либо порядку заполнения первого координатного угла в пространстве (x,y,t) ячейками, при котором новая ячейка прикладывается тремя гранями к ранее уложенным ячейкам или координатным плоскостям.


Пример 9.2. Построить разностные схемы для решения первой начально-краевой задачи для двумерного уравнения теплопроводности (где \Omega= (0,L)\times (0,M),~ \Gamma — граница области \Omega;~ \psi(x,y),~ \varphi(y,t) — заданные функции; L,\,M,\,T — заданные числа):

u_t=u_{xx}+u_{yy},~~ (x,y,t)\in (0,L)\times (0,M)\times (0,T),

u(x,y,0)= \psi(x,y),~~ (x,y)\in [0,L]\times [0,M], (начальное условие)

u(x,y,t)= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T, (краевое условие)


▼ Решение, формулы (9.10)-(9.13)

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


Шеститочечный шаблон

Для аппроксимации производных используем первую формулу в (9.5) и четвертые формулы в (9.6),(9.7):


\begin{gathered}\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n}}{\tau}= \frac{\widehat{u}_{i+1,j}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h_y^2}\,,\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1. \end{gathered}
(9.10)

Начальные и краевые условия представляются в виде


\begin{gathered}\widehat{u}_{i,j}^{0}=\psi_{i,j},\quad 0 \leqslant i \leqslant I,\quad 0 \leqslant j \leqslant J,\\[2pt] \widehat{u}_{0,j}^{n}= \varphi_{0,j}^{n},\quad \widehat{u}_{I,j}^{n}= \varphi_{I,j}^{n},\quad 0 \leqslant j \leqslant J,\quad 1 \leqslant n \leqslant N,\\[2pt] \widehat{u}_{i,0}^{n}= \varphi_{i,0}^{n},\quad \widehat{u}_{i,J}^{n}= \varphi_{i,J}^{n},\quad 0 \leqslant i \leqslant I,\quad 1 \leqslant n \leqslant N,\\[2pt] \end{gathered}
(9.11)

где \varphi_{i,j}^{n}= \varphi(x_{i},y_{j},t^{n}),~ \psi_{i,j}= \psi(x_{i},y_{j}). Начальные условия задают значения искомой функции на нижнем основании параллелепипеда, а краевые условия — на боковых гранях.


Соотношения (9.10), (9.11) образуют явную двухслойную разностную схему. Можно показать, что она условно устойчива (устойчива при выполнении условия \tau(h_{x}^{-2}+ h_{y}^{-2}) \leqslant 1\!\!\not{\phantom{|}}\,2). При h_x= h_y=h имеем условие \frac{\tau}{h^2} \leqslant \frac{1}{4}. В этом случае (9.10) принимает вид


\begin{gathered}\begin{aligned}\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n} }{\tau}&= \frac{\widehat{u}_{i+1,j}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h^2}+ \frac{\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h^2}=\\ &=\frac{\widehat{u}_{i+1,j}^{n}+ \widehat{u}_{i-1,j}^{n}+ \widehat{u}_{i,j-1}^{n}+ \widehat{u}_{i,j+1}^{n}-4\,\widehat{u}_{i,j}^{n}}{h^2}\,,\end{aligned}\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1.\end{gathered}

Из этого соотношения следует, что при фиксированных i,\,k,\,n искомое значение \widehat{u}_{i,j}^{n+1} определяется явным образом через известные значения на n-м слое:


\widehat{u}_{i,j}^{n+1}= \widehat{u}_{i,j}^{n}+ \frac{\tau}{h^2} \bigl[\widehat{u}_{i+1,j}^{n}+ \widehat{u}_{i-1,j}^{n}+ \widehat{u}_{i,j-1}^{n}+ \widehat{u}_{i,j+1}^{n}-4\,\widehat{u}_{i,j}^{n}\bigr].

Построим неявную схему, применяя шеститочечный шаблон, изображенный на рис. 9.4,б. Для аппроксимации производных используем первую формулу в (9.5) и четвертые в (9.6),(9.7), записанные на (n+1)-м временном слое.


\begin{gathered}\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n}}{\tau}= \frac{\widehat{u}_{i+1,j}^{n+1}-2\,\widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i-1,j}^{n+1}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n+1}-2\,\widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i,j-1}^{n+1} }{h_y^2}\,,\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1. \end{gathered}
(9.12)

Соотношения (9.12),(9.11) образуют неявную двухслойную схему. Можно показать, что она является абсолютно устойчивой. Соотношение (9.12), записанное для всех внутренних узлов (n+1)-го слоя, порождает систему линейных алгебраических уравнений для определения значений искомой функции в узлах. Каждое уравнение этой системы содержит только пять неизвестных. Следовательно, матрица системы является пятидиагональной и сильно разреженной. Для решения системы применяются, как правило, итерационные методы и метод прогонки.


Десятиточечный шаблон

Получим еще одну неявную схему, используя десятиточечный шаблон, изображенный на рис. 9.5. Он является обобщением шеститочечного шаблона (см. рис. 8.8,6) на двумерный случай. Аппроксимация уравнения имеет вид


\begin{gathered}\begin{aligned}\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n}}{\tau}&= \lambda\! \left(\frac{\widehat{u}_{i+1,j}^{n+1}-2\,\widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i-1,j}^{n+1}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n+1}-2\,\widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i,j-1}^{n+1} }{h_y^2}\right)\!,\\ &\quad +(1-\lambda)\! \left(\frac{\widehat{u}_{i+1,j}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n} }{h_y^2}\right)\!,\end{aligned}\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1, \end{gathered}
(9.13)

где 0 \leqslant \lambda \leqslant 1 — параметр схемы. При \lambda=0 схема становится явной и значение решения на (n+1)-м слое непосредственно вычисляется по известным значениям с предыдущего слоя, а при \lambda=1 совпадает с неявной схемой (9.12),(9.11). Начальные и краевые условия записываются в виде (9.11). Соотношения (9.13),(9.11) образуют неявную двухслойную разностную схему. Можно доказать, что схема имеет второй порядок аппроксимации по h_x,\,h_y, а по \tau — первый порядок при \lambda\ne \frac{1}{2} — и второй при \lambda= \frac{1}{2}. Схема является условно устойчивой (устойчивой при \lambda \geqslant \frac{1}{2}-\frac{1}{4\tau}\! \left(\frac{1}{h_x^2}+ \frac{1}{h_y^2}\right)^{-1}).


Пример 9.3. Построить разностные схемы для первой начально-краевой задачи для двумерного волнового уравнения (где \Omega= (0,L)\times (0,M),~ \Gamma — граница области \Omega;~ \psi_1(x,y),\,\psi_2(x,y),\, \varphi(x,y,t) — заданные функции):

u_{tt}= c^2(u_{xx}+u_{yy}),~~ (x,y,t)\in \Omega\times (0,T),

u(x,y,0)=\psi_1(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)

\frac{\partial u(x,y,0)}{\partial t}=\psi_2(x,y),~~ (x,y)\in \Omega\cup \Gamma, (начальное условие)

u(x,y,t)= \varphi(x,y,t),~~ (x,y)\in \Gamma,~ 0<t \leqslant T (краевое условие)


▼ Решение, формулы (9.14)-(9.15)

Сначала построим явную схему. Для этого воспользуемся шаблоном, изображенным на рис. 9.6,а, обобщающим шаблон «крест» , применяемый для решения одномерного волнового уравнения. Для аппроксимации частных производных будем использовать четвертые формулы в (9.5)–(9.7):


\begin{gathered}\frac{u_{i,j}^{n+1}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{\tau^2}= \frac{\widehat{u}_{i+1,j}^{n}-2\, \widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n}-2\, \widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h_y^2}\,,\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1. \end{gathered}
(9.14)

Функциональные начальное и краевое условия представляются так же, как и в примере 9.1. Дифференциальное начальное условие следует аппроксимировать со вторым порядком, чтобы не ухудшать свойств схемы (см. пример 8.6). В результате имеем явную трехслойную схему, имеющую второй порядок аппроксимации и условно устойчивую (устойчивую при выполнении условия \tau< \bigl(h_{x}^{-2}+ h_{y}^{-2}\bigr)^{-1\!\not{\phantom{|}}\,\,2}).


Пятнадцатиточечный шаблон

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


\begin{gathered}\begin{aligned}\frac{u_{i,j}^{n+1}-2u_{i,j}^{n}+u_{i,j}^{n-1}}{\tau^2}&= \lambda\! \left(\frac{\widehat{u}_{i+1,j}^{n+1}-2\, \widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i-1,j}^{n+1}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n+1}-2\, \widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i,j-1}^{n+1}}{h_y^2}\right)+\\ &\quad +(1-\lambda)\! \left(\frac{\widehat{u}_{i+1,j}^{n}-2\, \widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n}-2\, \widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h_y^2}\right)+\\ &\quad +\lambda\! \left(\frac{\widehat{u}_{i+1,j}^{n-1}-2\, \widehat{u}_{i,j}^{n-1}+ \widehat{u}_{i-1,j}^{n-1}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}^{n-1}-2\, \widehat{u}_{i,j}^{n-1}+ \widehat{u}_{i,j-1}^{n-1}}{h_y^2}\right)\!,\end{aligned}\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\quad 1 \leqslant n \leqslant N-1. \end{gathered}
(9.15)

Схема (9.15), (9.11) имеет второй порядок аппроксимации и является трехслойной. При \lambda=0 она переходит в явную. При \frac{1}{4} \leqslant \lambda \leqslant \frac{1}{2} схема абсолютно устойчива, но относится к неэкономичным.




Метод расщепления


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


Пример 9.4. Построить схему метода расщепления для решения задачи Коши (где \psi(x,y) — заданная функция):

u_{t}= u_{xx}+ u_{yy},~~ (x,y,t)\in \mathbb{R}^2\times (0,T),

u(x,y,0)= \psi(x,y),~~ (x,y)\in \mathbb{R}^2 (начальное условие).


Введем в рассмотрение дифференциальный оператор \mathcal{A}u \equiv u_{xx}+ u_{yy}. Тогда решаемое уравнение записывается в операторной форме: u_{t}= \mathcal{A}u. Предположим, что решение u(x,y,t^n) поставленной задачи известно в момент времени t=t^n. Выразим решение u(x,y,t^{n-1}) при t=t^{n+1} для фиксированных x и y через решение u(x,y,t^n), пользуясь формулой Тейлора:


\begin{aligned}u(x,y,t^{n+1})&= u(x,y,t^n+\tau)= u(x,y,t^n)+ \frac{\tau}{1!}\cdot \frac{\partial u(x,y,t^n)}{\partial t}+ O(\tau^2)=\\ &=u(x,y,t^n)+ \frac{\tau}{1!}\cdot \mathcal{A} u(x,y,t^n)+ O(\tau^2)= (\mathcal{E}+ \tau \mathcal{A}) u(x,y,t^n)+ O(\tau^2), \end{aligned}

где производная \frac{\partial u(x,y,t^n)}{\partial t} во втором слагаемом заменена правой частью исходного дифференциального уравнения, записанного в операторной форме, а \mathcal{E} — тождественный оператор \bigl(\mathcal{E}u(x,y,t^n)= u(x,y,t^n)\bigr).


Зададим сетку D_h= \{x_{i}, y_{j}, t^n\}, где x_{i}= ih,~ y_{j}= jh,~ t^n=n\tau,~ i,j=0,\pm1,\pm2,\ldots;~ n=0,1,\ldots,N;~ T=N\tau. Оператор правой части исходного уравнения представим в виде суммы двух операторов:


\mathcal{A}= \mathcal{A}_1+ \mathcal{A}_2,\quad \mathcal{A}_1 \equiv u_{xx},\quad \mathcal{A}_2 \equiv u_{yy}.

Наряду с исходной задачей рассмотрим две вспомогательные задачи Коши:


– задача 1, в которой решается уравнение


\frac{\partial \nu(x,y,t)}{\partial t}= \frac{\partial^2 \nu (x,y,t)}{\partial x^2}\,,\quad (x,y,t)\in \mathbb{R}^2\times (t^n,t^{n+1}),

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

\begin{aligned}& \nu(x,y,t^n)= u(x,y,t^n), & & (x,y)\in \mathbb{R}^2, & & n=1,2,\ldots,N-1,\\ & \nu(x,y,t^0)= \psi(x,y), & & (x,y)\in \mathbb{R}^2, & & n=0; \end{aligned}

– задача 2, в которой решается уравнение

\frac{\partial w(x,y,t)}{\partial t}= \frac{\partial^2 w(x,y,t)}{\partial y^2}\,,\quad (x,y,t)\in \mathbb{R}^2\times (t^n,t^{n+1}),

с начальным условием w(x,y,t^{n})= \nu(x,y,t^{n+1}),~ (x,y)\in \mathbb{R}^2,~ n=0,1,\ldots, N-1.


Поставленные задачи могут быть решены последовательно: сначала задача 1, а затем задача 2. Решая задачу 1, можно найти функцию \nu(x,y,t^{n+1}) на (n+1)-м слое, которая определяет начальное условие задачи 2. Решая задачу 2, можно найти функцию w(x,y,t^{n+1}) на том же слое. Задача 1 при фиксированном y и задача 2 при фиксированном x фактически являются одномерными. Установим связь между функциями \nu(x,y,t^{n+1}),~ w(x,y,t^{n+1}) и решением u(x,y,t^{n+1}) исходной задачи на (n+1)-м слое. Пользуясь формулой Тейлора, для функции \nu(x,y,t^{n+1}) получим соотношение \nu(x,y,t^{n+1})= (\mathcal{E}+ \tau \mathcal{A}_1) \nu(x,y,t^{n})+ O(\tau^2)= (\mathcal{E}+ \tau \mathcal{A}_1) u(x,y,t^{n})+ O(\tau^2), где функция \nu(x,y,t^{n}) заменена на u(x,y,t^{n}) в силу начального условия в задаче 1. Для функции w(x,y,t^{n+1}) можно получить аналогичное выражение


w(x,y,t^{n+1})= (\mathcal{E}+ \tau \mathcal{A}_2) w(x,y,t^{n})+ O(\tau^2)= (\mathcal{E}+ \tau \mathcal{A}_2) \nu(x,y,t^{n})+ O(\tau^2);

\begin{aligned}w(x,y,t^{n+1})&= (\mathcal{E}+ \tau \mathcal{A}_2) \bigl[(\mathcal{E}+ \tau \mathcal{A}_1) u(x,y,t^{n})+ O(\tau^2)\bigr]+ O(\tau^2)=\\[4pt] &=\mathcal{E}^2 u(x,y,t^{n})+ \tau(\mathcal{A}_1+ \mathcal{A}_2)u(x,y,t^{n})+ \tau^2 \mathcal{A}_2 \mathcal{A}_1 u(x,y,t^{n})+ O(\tau^2)=\\[4pt] &=(\mathcal{E}+ \tau \mathcal{A}) u(x,y,t^{n})+ O(\tau^2)= u(x,y,t^{n+1})+ O(\tau^2), \end{aligned}
(9.16)

так как \mathcal{E}^2 u(x,y,t^{n})= \mathcal{E} u(x,y,t^{n}),~ \tau^2 \mathcal{A}_2 \mathcal{A}_1 u(x,y,t^{n})= O(\tau^2). Таким образом, в результате последовательного решения задач 1 и 2 при t=t^{n+1} получаем функцию w(x,y,t^{n+1}), которая отличается от решения исходной задачи u(x,y,t^{n+1}) лишь на величину O(\tau^2). Поэтому можно положить u(x,y,t^{n+1})\cong w(x,y,t^{n+1}).


Следовательно, здесь процесс решения исходной двумерной задачи заменен процессом решения двух задач в силу расщепления дифференциального оператора \mathcal{A} на сумму двух операторов \mathcal{A}_1+ \mathcal{A}_2.


Замечания


1. В описываемом методе расщепление произведено на уровне задач. В методе переменных направлений расщепление на уровне задач не делается, но слагаемые оператора \mathcal{A} аппроксимируются на различных временных слоях.

2. Имеется модификация данного метода, в которой расщепление осуществляется по физическим процессам, когда отдельные слагаемые оператора с/1 соответствуют различным физическим явлениям (например, неравновесным химическим или колебательным процессам, процессам подогрева, двухфазным процессам и т.д.).


Многоточечные шаблоны

Сформируем две схемы расщепления, заменяя задачи 1 и 2 разностными аппроксимациями. При построении явной схемы применим шаблоны, изображенные на рис. 9.7. Для аппроксимации частных производных по времени применим первую формулу в (9.5), а для аппроксимации пространственных производных — последние формулы в (9.6) и (9.7), имеющие второй порядок аппроксимации. Тогда для задачи 1 получаем (шаблон на рис. 9.7,а), где \psi_{i,j}= \psi(x_{i},y_{j})\colon


\begin{gathered}\frac{\widehat{\nu}_{i,j}^{n+1}-\widehat{\nu}_{i,j}^{n}}{\tau}= \frac{\widehat{\nu}_{i+1,j}^{n}-2\,\widehat{\nu}_{i,j}^{n}+ \widehat{\nu}_{i-1,j}^{n}}{h^2}\,,\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1;\\[4pt] \widehat{\nu}_{i,j}^{0}= \psi_{i,j},\quad i,j=0,\pm1,\pm2,\ldots;\\[4pt] \widehat{\nu}_{i,j}^{n}= \widehat{u}_{i,j}^{n},\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1. \end{gathered}
(9.17)

Для задачи 2 аналогично получаем (шаблон на рис. 9.7,б)


\begin{gathered}\frac{\widehat{w}_{i,j}^{n+1}-\widehat{w}_{i,j}^{n}}{\tau}= \frac{\widehat{w}_{i,j+1}^{n}-2\,\widehat{w}_{i,j}^{n}+ \widehat{w}_{i,j-1}^{n}}{h^2}\,,\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1;\\[4pt] \widehat{w}_{i,j}^{n}= \widehat{\nu}_{i,j}^{n+1},\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1. \end{gathered}
(9.18)

Заметим, что (9.17) с учетом условия \widehat{\nu}_{i,j}^{n}= \widehat{u}_{i,j}^{n} можно переписать в форме


\widehat{\nu}_{i,j}^{n+1}= \widehat{u}_{i,j}^{n}+ \tau\, \frac{\widehat{u}_{i+1,j}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i-1,j}^{n}}{h^2}, а (9.18) — в виде \widehat{w}_{i,j}^{n+1}= \widehat{\nu}_{i,j}^{n+1}+ \tau\, \frac{\widehat{w}_{i,j+1}^{n}-2\,\widehat{w}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h^2}\,.

С учетом (9.16) и того, что начальные условия при n=0 аппроксимируются точно, эта явная двухслойная разностная схема аппроксимирует исходную задачу с порядком O(\tau^2+ h^2). Схема является условно устойчивой (устойчивой при \tau \leqslant h^2\!\!\not{\phantom{|}}\,2).




Методика вычислений по явной схеме


1. Задать значения шагов сетки h,\,\tau так, что N\,\tau=T, где N — целое положительное число, причем значение шага по времени выбрать из условия устойчивости \tau \leqslant h^2\!\!\not{\phantom{|}}\,2. Вычислить \widehat{\nu}_{i,j}^{0}= \widehat{u}_{i,j}^{0}= \psi_{i,j}, i,j= 0,\pm1, \pm2, \ldots, где \psi_{i,j}= \psi(x_{i},y_{j}). Положить n=0.


2. Найти решение на (n+1)-м слое.


Первый этап. Решить задачу 1 при каждом j=0,\pm1, \pm2, \ldots\colon


\widehat{\nu}_{i,j}^{n+1}= \widehat{u}_{i,j}^{n}= \tau\cdot \frac{\widehat{u}_{i+1,j}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{\nu}_{i-1,j}^{n}}{h^2}\,,\quad i= 0,\pm1, \pm2, \ldots

Второй этап. Решить задачу 2 при каждом i=0,\pm1, \pm2, \ldots\colon


\widehat{w}_{i,j}^{n+1}= \widehat{\nu}_{i,j}^{n+1}= \tau\cdot \frac{\widehat{w}_{i,j+1}^{n}-2\,\widehat{w}_{i,j}^{n}+ \widehat{w}_{i,j-1}^{n}}{h^2}\,,\quad j= 0,\pm1, \pm2, \ldots

Положить \widehat{u}_{i,j}^{n+1}= \widehat{w}_{i,j}^{n+1},~ i,j= 0,\pm1, \pm2, \ldots


3. Если n=N, вычисления завершить. Иначе положить n=n+1 и перейти к пункту 2.


При построении неявной схемы применим шаблоны, изображенные на рис. 9.8. Для аппроксимации частных производных по времени применим первую формулу в (9.5), а для аппроксимации пространственных производных — последние формулы в (9.6) и (9.7), записанные на (n+1)-м временном слое.


Шеститочечные шаблоны

Тогда для задачи 1 получаем (шаблон на рис. 9.8,а), где \psi_{i,j}= \psi(x_{i},y_{j})\colon


\begin{gathered}\frac{\widehat{\nu}_{i,j}^{n+1}-\widehat{\nu}_{i,j}^{n}}{\tau}= \frac{\widehat{\nu}_{i+1,j}^{n+1}-2\,\widehat{\nu}_{i,j}^{n+1}+ \widehat{\nu}_{i-1,j}^{n+1}}{h^2}\,,\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1;\\[4pt] \widehat{\nu}_{i,j}^{0}= \psi_{i,j},\quad i,j=0,\pm1,\pm2,\ldots;\\[4pt] \widehat{\nu}_{i,j}^{n}= \widehat{u}_{i,j}^{n},\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1. \end{gathered}
(9.19)

Для задачи 2 аналогично получаем (шаблон на рис. 9.8,б)


\begin{gathered}\frac{\widehat{w}_{i,j}^{n+1}-\widehat{w}_{i,j}^{n}}{\tau}= \frac{\widehat{w}_{i,j+1}^{n+1}-2\,\widehat{w}_{i,j}^{n+1}+ \widehat{w}_{i,j-1}^{n+1}}{h^2}\,,\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1;\\[4pt] \widehat{w}_{i,j}^{n}= \widehat{\nu}_{i,j}^{n+1},\quad i,j=0,\pm1,\pm2,\ldots;\quad n=0,1,\ldots,N-1. \end{gathered}
(9.20)

Приведенные соотношения образуют абсолютно устойчивую двухслойную неявную разностную схему метода расщепления, имеющую порядок аппроксимации O(\tau^2+h^2).




Методика вычислений по неявной схеме


1. Задать значения шагов сетки h,\,\tau так, что N\,\tau=T, где N — целое положительное число, причем значения шагов выбрать из условия достижения заданной точности. Вычислить \widehat{\nu}_{i,j}^{0}= \widehat{u}_{i,j}^{0}= \psi_{i,j}, где \psi_{i,j}= \psi(x_{i},y_{j}). Положить n=0.


2. Найти решение на (n+1)-м слое.


Первый этап. Решить задачу 1 при каждом j=0,\pm1, \pm2, \ldots \colon


\frac{1}{h^2}\cdot \widehat{\nu}_{i+1,j}^{n+1}-\left(\frac{2}{h^2}+ \frac{1}{\tau}\right)\!\cdot \widehat{\nu}_{i,j}^{n+1}+ \frac{1}{h^2}\cdot \widehat{\nu}_{i-1,j}^{n+1}=-\frac{1}{\tau}cdot \widehat{u}_{i,j}^{n},\quad i=0,\pm1,\pm2,\ldots

применяя метод прогонки для системы линейных алгебраических уравнений с трехдиагональной матрицей (система получена из (9.19) при \widehat{\nu}_{i,j}^{n}= \widehat{u}_{i,j}^{n}).


Второй этап. Решить задачу 2 при каждом j=0,\pm1,\pm2, \ldots\colon


\frac{1}{h^2}\cdot \widehat{w}_{i,j+1}^{n+1}-\left(\frac{2}{h^2}+ \frac{1}{\tau}\right)\!\cdot \widehat{w}_{i,j}^{n+1}+ \frac{1}{h^2}\cdot \widehat{w}_{i,j-1}^{n+1}=-\frac{1}{\tau}cdot \widehat{\nu}_{i,j}^{n+1},\quad j=0,\pm1,\pm2,\ldots

применяя метод прогонки для системы линейных алгебраических уравнений с трехдиагональной матрицей (система получена из (9.20) при \widehat{w}_{i,j}^{n}= \widehat{\nu}_{i,j}^{n+ 1}). Положить \widehat{u}_{i,j}^{n+1}= \widehat{w}_{i,j}^{n+1},~ i,j= 0,\pm1, \pm2, \ldots.


3. Если n=N, вычисления завершить. Иначе положить n=n+1 и перейти к пункту 2.


Замечания

1. Для решаемых на первом и втором этапах методики систем линейных алгебраических уравнений выполняется условие преобладания диагональных элементов, следовательно, метод прогонки устойчив.

2. Для пояснения общей схемы вычислений здесь рассматривалась задача Коши, поэтому индексы i и j были не ограничены. Однако эти схемы могут быть легко приспособлены для решения смешанных краевых задач с ограниченным множеством D.




Метод переменных направлений


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


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

u_{t}= u_{xx}+ u_{yy}+ f(x,y,t),\quad (x,y,t)\in (0,L)\times (0,M)\times (0,T),

u(x,y,0)= \psi(x,y,t),\quad (x,y)\in [0,L]\times [0,M], (начальное условие)

u(x,y,t)= \varphi(x,y,t),\quad (x,y)\in \Gamma,~ 0<t \leqslant T (краевое условие).


где \Omega= (0,L)\times (0,M),~ \Gamma — граница области \Omega (прямоугольника); f(x,y,t) — плотность источников тепла; \psi(x,y),\, \varphi(x,y,t) — заданные функции; L,\,M,\,T — заданные числа.


Как и в методе расщепления, исходное уравнение аппроксимируется совокупностью двух разностных схем, каждая из которых соответствует только одному пространственному направлению. Однако оператор \mathcal{A} здесь не расщепляется на сумму \mathcal{A}_1+ \mathcal{A}_2 двух операторов на уровне задач, а каждый элемент суммы аппроксимируется с помощью явных и неявных конструкций. Для этого наряду со слоем t=t^{n+1}, расчет которого осуществляется на данном этапе, вводится один дополнительный (промежуточный) слой t=t^{n+1\!\not{\phantom{|}}\,\,2}. Тогда для перехода от слоя t=t^n к слою t=t^{n+1} с использованием слоя t=t^{n+ 1\!\not{\phantom{|}}\,\,2} исходное дифференциальное уравнение аппроксимируется двумя разностными уравнениями, одно из которых связывает слои t=t^{n} и t=t^{n+ 1\!\not{\phantom{|}}\,\,2}, а второе — слои t=t^{n+ 1\!\not{\phantom{|}}\,\,2} и t=t^{n+1}.


При реализации метода используем девятиточечный шаблон, изображенный на рис. 9.9. Переход с n-го временного слоя на (n+1)-й разбивается на два этапа.


Девятиточечный шаблон

Первый этап. Переход на дополнительный промежуточный (n+1/2)-й слои с шагом \tau\!\!\not{\phantom{|}}\,2. Частная производная по времени аппроксимируется по первой формуле в (9.5), производная по x аппроксимируется на (n+1/2)-м слое, производная по y — на n-м слое. Входящая в правую часть исходного уравнения функция f(x,y,t) заменяется своим сеточным представлением. Соответствующая разностная схема имеет вид:


\begin{gathered}\frac{\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-\widehat{u}_{i,j}^{n} }{\tau\!\!\not{\phantom{|}}\,2}= \frac{\widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-2\, \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}+ \widehat{u}_{i-1,j}^{n+ 1\!\not{\phantom{|}}\,\,2} }{h_{x}^{2}}+ \frac{\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}}{h_{y}^{2}}+ f_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\end{gathered}
(9.21)

где f_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}= f\! \left(x_{i},y_{j}, t^n+ \frac{1}{2}\tau\right), или в традиционной записи:


\frac{\tau}{2h_{x}^{2}}\,\widehat{u}_{i-1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-\left(1+ \frac{\tau}{h_{x}^{2}}\right)\! \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}+ \frac{\tau}{2h_{x}^{2}}\, \widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}= F_{i,j}^{n},\quad 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,

где F_{i,j}^{n}=-\frac{\tau}{2}f_{i.j}^{n+1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{2h_{y}^{2}}\bigl(\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}\bigr)-\widehat{u}_{i,j}^{n}. Заметим, что вместо значения f_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2} в (9.21) можно использовать f_{i,j}^{n}.


Для каждого j=1,\ldots,J-1 требуется решить трехдиагональную систему линейных алгебраических уравнений, так как каждое уравнение содержит три неизвестных значения \widehat{u}_{i-1,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\, \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\, \widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}, а остальные значения берутся с n-го слоя или вычисляются, например f_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}. Иными словами, при таком переходе схема является неявной по направлению x и явной по направлению y. Искомые значения на промежуточном слое вычисляются методом прогонки по направлению x, т.е. продольному направлению.


Второй этап. Переход на (n+1)-й временной слой с промежуточного (n+1/2)-го с шагом \tau/2. Частная производная по времени аппроксимируется по первой формуле в (9.5), производная по x аппроксимируется на (n+1/2)-м слое, производная по y — на (n+1)-м слое. Соответствующая разностная схема имеет вид:


\begin{gathered}\frac{\widehat{u}_{i,j}^{n+1}-\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2} }{\tau\!\!\not{\phantom{|}}\,2}= \frac{\widehat{u}_{i+ 1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-2\, \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}+ \widehat{u}_{i-1,j}^{n+ 1\!\not{\phantom{|}}\,\,2} }{h_{x}^{2}}+ \frac{\widehat{u}_{i,j+1}^{n+1}-2\,\widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i,j-1}^{n+1}}{h_{y}^{2}}+ f_{i,j}^{n+1},\\[2pt] 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,\end{gathered}
(9.22)

или в традиционной записи:


\frac{\tau}{2h_{y}^{2}}\,\widehat{u}_{i,j-1}^{n+ 1}-\left(1+ \frac{\tau}{h_{y}^{2}}\right)\! \widehat{u}_{i,j}^{n+1}+ \frac{\tau}{2h_{y}^{2}}\, \widehat{u}_{i,j+1}^{n+1}= \Phi_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\quad 1 \leqslant i \leqslant I-1,\quad 1 \leqslant j \leqslant J-1,

где \Phi_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}=-\frac{\tau}{2} f_{i.j}^{n}-\frac{\tau}{2h_{x}^{2}}\bigl(\widehat{u}_{i-1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-2\,\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}+ \widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}\bigr)-\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}. Заметим, что вместо значения f_{i.j}^{n} в (9.22) можно использовать f_{i.j}^{n+ 1\!\not{\phantom{|}}\,\,2}.


Для каждого значения i=1,2,\ldots,I-1 требуется решить трехдиагональную систему линейных алгебраических уравнений, так как каждое уравнение содержит три неизвестных значения \widehat{u}_{i,j-1}^{n+1},\, \widehat{u}_{i,j}^{n+1},\, \widehat{u}_{i,j+ 1}^{n+1}, а остальные значения берутся с (n+1/2)-го слоя. Иными словами, при таком переходе схема является явной по направлению x и неявной по направлению y. Искомые значения на целом слое вычисляются методом прогонки по направлению y, т.е. в поперечном направлении. Диагональные элементы в системах (9.21),(9.22) преобладают, следовательно, процесс прогонки устойчив, а решение разностной задачи существует и единственно. К достоинствам метода, достигнутым благодаря введению промежуточного слоя, следует отнести расщепление исходной задачи на две более простые, решаемые с помощью алгоритма скалярной прогонки.


Начальные и краевые условия представляются в виде (где \varphi_{i,j}^{n}= \varphi(x_{i}, y_{j}, t^{n}),~ \psi_{i,j}= \psi(x_{i},y_{j}))


\begin{aligned}& \widehat{u}_{i,j}^{0}= \psi_{i,j},\quad 0 \leqslant i \leqslant I,~ 0 \leqslant j \leqslant J,\quad \text{(nulevoy sloy)}\\[4pt] & \widehat{u}_{0,j}^{n}= \varphi_{0,j}^{n},~ \widehat{u}_{I,j}^{n}= \varphi_{I,j}^{n},\quad 0 \leqslant j \leqslant J,~ 1 \leqslant n \leqslant N,\quad \text{(na ploskostyah, perpendikulyarnih osi)}~Ox \\[4pt] & \widehat{u}_{i,0}^{n}= \varphi_{i,0}^{n},~ \widehat{u}_{i,J}^{n}= \varphi_{I,j}^{n},\quad 1 \leqslant i \leqslant I,~ 1 \leqslant n \leqslant N,\text{(na ploskostyah, perpendikulyarnih osi)}~Oy \end{aligned}
(9.23)

Для промежуточного слоя требуются значения \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2} на сторонах расчетной области, определяемых уравнениями x=0 и x=L. Вычитая (9.22) из (9.21), получаем


\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}= \frac{1}{2}\bigl(\widehat{u}_{i,j}^{n+ 1}+ \widehat{u}_{i,j}^{n}\bigr)-\frac{\tau}{4h_y^2} \bigl(\widehat{u}_{i,j+1}^{n+1}-2\, \widehat{u}_{i,j}^{n+1}+ \widehat{u}_{i,j+1}^{n+1}\bigr)+ \frac{\tau}{4hy^2} \bigl(\widehat{u}_{i,j+ 1}^{n}-2\, \widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j+1}^{n}\bigr).

Отсюда при i=0~(x=0) имеем


\widehat{u}_{0,j}^{n+ 1\!\not{\phantom{|}}\,\,2}= \frac{1}{2}\bigl(\widehat{u}_{0,j}^{n+ 1}+ \widehat{u}_{0,j}^{n}\bigr)-\frac{\tau}{4h_y^2} \bigl(\widehat{u}_{0,j+1}^{n+1}-2\, \widehat{u}_{0,j}^{n+1}+ \widehat{u}_{0,j+1}^{n+1}\bigr)+ \frac{\tau}{4hy^2} \bigl(\widehat{u}_{0,j+ 1}^{n}-2\, \widehat{u}_{0,j}^{n}+ \widehat{u}_{0,j+1}^{n}\bigr).

Такой способ обеспечивает второй порядок аппроксимации. Аналогично получается соотношение при i=I~(x=L). Схема (9.21)–(9.23) метода переменных направлений имеет второй порядок аппроксимации и является абсолютно устойчивой.


Методика расчетов


1. Задать сетку D_h= \bigl\{x_{i},y_{j},t^{n}\bigr\}, где x_{i}=i\,h_{x},\, 0 \leqslant i \leqslant I;~ y_{j}=j\,h_y,\, 0 \leqslant j \leqslant J;~ t^n=n\,\tau,\, 0 \leqslant n \leqslant N; N,\,J,\,I — целые положительные числа; h_x,\,h_y,\,\tau — величины шагов по пространственным координатам и времени (для простоты принимаются постоянными); I\,h_x=L,\, J\,h_y=M,\, N\,\tau=T.


2. Вычислить


\begin{gathered}\widehat{u}_{i,j}^{0}= \psi_{i,j},\quad 0 \leqslant i \leqslant I,\quad 0 \leqslant j \leqslant J,\\[2pt] \widehat{u}_{0,j}^{n}= \varphi_{0,j}^{n},\quad \widehat{u}_{I,j}^{n}= \varphi_{I,j}^{n},\quad 0 \leqslant j \leqslant J,\quad 1 \leqslant n \leqslant N,\\[2pt] \widehat{u}_{i,0}^{n}= \varphi_{i,0}^{n},\quad \widehat{u}_{i,J}^{n}= \varphi_{i,J}^{n},\quad 0 \leqslant i \leqslant I,\quad 1 \leqslant n \leqslant N. \end{gathered}

где \varphi_{i,j}^{n}= \varphi(x_{i},y_{j},t^{n}),~ \psi_{i,j}= \psi(x_{i},y_{j}). Положить n=0.


3. Найти решение на (n+1)-м слое.


Первый этап. Перейти на дополнительный промежуточный (n+1/2)-й слой с шагом \tau\!\!\not{\phantom{|}}\,2\colon для каждого j=1,\ldots,J-1 решить трехдиагональную систему линейных алгебраических уравнений


\frac{\tau}{2h_x^2}\cdot \widehat{u}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2}-\left(1+ \frac{\tau}{h_x^2}\right)\!\cdot \widehat{u}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}+ \frac{\tau}{2h_x^2}\cdot \widehat{u}_{i+1,j}^{n+1\!\not{\phantom{|}}\,\,2}= F_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},\quad 1 \leqslant i \leqslant I-1,

где F_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}=-\frac{\tau}{2}f_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{2h_y^2}\bigl(\widehat{u}_{i,j+1}^{n}-2\,\widehat{u}_{i,j}^{n}+ \widehat{u}_{i,j-1}^{n}\bigr)-\widehat{u}_{i,j}^{n},
\widehat{u}_{0,j}^{n+1\!\not{\phantom{|}}\,\,2}= \frac{1}{2}\bigl(\widehat{u}_{0,j}^{n+1}-\widehat{u}_{0,j}^{n}\bigr)-\frac{\tau}{4h_y^2} \bigl(\widehat{u}_{0,j+1}^{n+1}-2\,\widehat{u}_{0,j}^{n+1}+ \widehat{u}_{0,j+1}^{n+1}\bigr)+ \frac{\tau}{4h_y^2} \bigl(\widehat{u}_{0,j+1}^{n}-2\,\widehat{u}_{0,j}^{n}+ \widehat{u}_{0,j+1}^{n}\bigr),
\widehat{u}_{I,j}^{n+1\!\not{\phantom{|}}\,\,2}= \frac{1}{2}\bigl(\widehat{u}_{I,j}^{n+1}-\widehat{u}_{I,j}^{n}\bigr)-\frac{\tau}{4h_y^2} \bigl(\widehat{u}_{I,j+1}^{n+1}-2\,\widehat{u}_{I,j}^{n+1}+ \widehat{u}_{I,j+1}^{n+1}\bigr)+ \frac{\tau}{4h_y^2} \bigl(\widehat{u}_{I,j+1}^{n}-2\,\widehat{u}_{I,j}^{n}+ \widehat{u}_{I,j+1}^{n}\bigr).

Искомые значения \widehat{u}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2},\, \widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\, \widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2} на промежуточном слое вычисляются методом прогонки по направлению x.


Второй этап. Перейти на (n+1) -й временной слой с промежуточного (n+1/2)-го с шагом \tau\!\!\not{\phantom{|}}\,2\colon для каждого значения i=1,\ldots,I-1 решить трехдиагональную систему линейных алгебраических уравнений


\frac{\tau}{2h_y^2}\cdot \widehat{u}_{i,j-1}^{n+1\!\not{\phantom{|}}\,\,2}-\left(1+ \frac{\tau}{h_y^2}\right)\!\cdot \widehat{u}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}+ \frac{\tau}{2h_y^2}\cdot \widehat{u}_{i,j+1}^{n+1\!\not{\phantom{|}}\,\,2}= \Phi_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},\quad 1 \leqslant j \leqslant J-1,

где \Phi_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}=-\frac{\tau}{2}f_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{2h_x^2}\bigl(\widehat{u}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2}-2\,\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}+ \widehat{u}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}\bigr)-\widehat{u}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}. Искомые значения \widehat{u}_{i,j-1}^{n+1},\, \widehat{u}_{i,j}^{n+1},\, \widehat{u}_{i,j+ 1}^{n+ 1} на целом слое вычисляются методом прогонки по направлению y.


4. Если n=N, вычисления завершить. Иначе положить n=n+1 и перейти к пункту 3.


Замечания

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

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


Рассмотрим задачу Дирихле для уравнения Лапласа с двумя независимыми переменными (где \varphi(x,y)):


u_{xx}+u_{yy}=0,\quad (x,y)\in \Omega,

\Bigl.{u}\Bigr|_{\Gamma}= \varphi(x,y),\quad (x,y)\in \Gamma (краевое условие).

Наряду с этой стационарной задачей рассмотрим нестационарную задачу для двумерного уравнения теплопроводности с теми же краевыми условиями и произвольно заданными начальными условиями:


u_{t}=u_{xx}+u_{yy},\quad (x,y,t)\in \Omega\times (0,+\infty),
u(x,y,0)= \psi(x,y),\quad (x,y)\in \Omega\cup \Gamma (начальное условие),
\Bigl.{u(x,y,t)}\Bigr|_{\Gamma}= \varphi(x,y),\quad (x,y)\in \Gamma,~ 0<t<+\infty (краевое условие).

Таким образом, вводится дополнительная независимая переменная t и исходная задача усложняется, так как увеличивается ее размерность. Однако при этом становится возможным выбор наиболее экономичной разностной схемы решения уравнения теплопроводности. Можно доказать, что решение нестационарной задачи при t\to+\infty сходится к решению стационарной задачи, а при достаточно большом t=T может быть принято за приближенное решение задачи Дирихле. Такой метод называется методом установления. Для прекращения вычислений можно использовать, например, условие \max_{i,j} \bigl| \widehat{u}_{i,j}^{n+l}-\widehat{u}_{i,j}^{n}\bigr| \leqslant \varepsilon, где \varepsilon> 0,~ l — положительное целое число (порядка нескольких десятков).


Для решения нестационарной задачи можно использовать схемы (9.21),(9.22) и метод переменных направлений. В последнем случае существует наилучший шаг по времени: \tau\cong \frac{1}{\pi} \sqrt{h_x^2+ h_y^2}{\left(\frac{1}{L^2}+ \frac{1}{M^2}\right)\!}^{-1\!\not{\phantom{|}}\,\,2}, при котором нестационарное решение выйдет на стационарное за наименьшее число шагов.




Метод дробных шагов


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


Пример 9.6. Построить схему метода дробных шагов для решения задачи Коши (где \psi(x,y) — заданная функция):
u_{t}=u_{xx}+u_{yy},\quad (x,y,t)\in \mathbb{R}^2\times (0,T),
u(x,y,0)=\psi(x,y),\quad (x,y)\in \mathbb{R} (начальное условие).

Зададим сетку D_h=\{x_{i},y_{j},t^n\}, где x_{i}=ih,~ y_{j}=jh,~ t^n= n\tau,~ i,j=0,\pm1,\pm2,\ldots,~ n=0,1,\ldots,N;~ T=N\tau.


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


\mathcal{A}= \mathcal{A}_1+\mathcal{A}_2,\quad \mathcal{A}_1 \equiv u_{xx},\quad \mathcal{A}_2 \equiv u_{yy}.

Далее они аппроксимируются отдельно друг от друга. В отличие от метода расщепления вводится вспомогательный (промежуточный) слой, соответствующий значению t=t^{n+ 1\!\not{\phantom{|}}\,\,2}. Как в методе переменных направлений, сначала осуществляется переход на промежуточный (n+1/2)-й слой с дробным шагом \tau\!\!\not{\phantom{|}}\,2, а затем с тем же шагом на целый (n+1)-й слой. Наряду с исходной задачей рассмотрим две вспомогательные задачи Коши:


задача 1, в которой решается уравнение

\frac{\partial \nu(x,y,t)}{\partial t}= \frac{\partial^2 \nu(x,y,t)}{\partial x^2},\quad (x,y,t)\in \mathbb{R}^2\times (t^n, t^{n+ 1\!\not{\phantom{|}}\,\,2}),

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

\begin{aligned}&\nu(x,y,t^n)= u(x,y,t^n), & & (x,y)\in \mathbb{R}^2, & & n=1,2,\ldots,N-1,\\ &\nu(x,y,t^0)= \psi(x,y), & & (x,y)\in \mathbb{R}^2, & & n=0;\end{aligned}

задана 2, в которой решается уравнение

\frac{\partial w(x,y,t)}{\partial t}= \frac{\partial^2 w(x,y,t)}{\partial x^2},\quad (x,y,t)\in \mathbb{R}^2\times (t^{n+ 1\!\not{\phantom{|}}\,\,2},t^{n+1}),

с начальным условием w(x,y,t^{n+ 1\!\not{\phantom{|}}\,\,2})= \nu(x,y,t^{n+ 1\!\not{\phantom{|}}\,\,2}),~ (x,y)\in \mathbb{R}^2,~ n=0,1,\ldots,N-1.


Поставленные задачи могут быть решены последовательно: сначала задача 1, а затем задача 2. Решая задачу 1, можно найти функцию \nu(x,y,t^{n+ 1\!\not{\phantom{|}}\,\,2}), которая определяет начальное условие задачи 2. Решая задачу 2, можно найти функцию w(x,y, t^{n+1}). Задача 1 при фиксированном y и задача 2 при фиксированном x фактически являются одномерными. Как и в методе расщепления, функция w(x,y,t^{n+1}) является приближенным решением исходной задачи на (n+1)-м слое, т.е. можно положить u(x,y,t^{n+1})\cong w(x,y,t^{n+1}). При реализации метода используем семиточечный шаблон, изображенный на рис. 9.10. Переход с n-го временного слоя на (n+1)-й разбивается на два этапа.


Семиточечный шаблон

Первый этап. Переход на дополнительный промежуточный (n+1/2)-й слой с шагом \tau\!\!\not{\phantom{|}}\,2. В задаче 1 частная производная по времени аппроксимируется по первой формуле в (9.5), производная по x аппроксимируется по последней формуле в (9.6) на (n+1/2)-м слое. Соответствующая разностная схема имеет вид:


\begin{gathered}\frac{\widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}-\widehat{\nu}_{i,j}^{n} }{\tau\!\!\not{\phantom{|}}\,2}= \frac{\widehat{\nu}_{i+1,j}^{n+ 1\!\not{\phantom{|}}\,\,2}-2\, \widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}+ \widehat{\nu}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2} }{h_x^2}\,,\quad i,j=0,\pm1,\pm2,\ldots,\quad n=0,1,\ldots,N-1;\\[2pt] \widehat{\nu}_{i,j}^{0}= \psi_{i,j},\quad i,j=0,\pm1,\pm2,\ldots,\\[2pt] \widehat{\nu}_{i,j}^{n}= \widehat{u}_{i,j}^{n},\quad i,j=0,\pm1,\pm2,\ldots,\quad n=1,2,\ldots,N-1, \end{gathered}
(9.24)

где \psi_{i,j}= \psi(x_{i},y_{j}), или в традиционной записи:


\frac{\tau}{2h_x^2}\,\widehat{\nu}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2}-\left(1+ \frac{\tau}{h_x^2}\right)\!\widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}+ \frac{\tau}{2h_x^2}\, \widehat{\nu}_{i+1,j}^{n+1\!\not{\phantom{|}}\,\,2}=-\widehat{u}_{i,j}^{n},\quad i,j=0,\pm1,\pm2,\ldots

Для каждого j требуется решить трехдиагональную систему линейных алгебраических уравнений, так как каждое уравнение содержит три неизвестных значения \widehat{\nu}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2},\, \widehat{\nu}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\, \widehat{\nu}_{i+1, j}^{n+ 1\!\not{\phantom{|}}\,\,2}, а остальные значения берутся с n-го слоя. Иными словами, при таком переходе схема является неявной. Искомые значения на промежуточном слое вычисляются методом прогонки по направлению x.


Второй этап. Переход на (n+1)-й временной слой с промежуточного (n+1/2)-го с шагом \tau\!\!\not{\phantom{|}}\,2. Частная производная по времени аппроксимируется по первой формуле в (9.5), производная по у аппроксимируется по последней формуле в (9.7) на (n+1)-м слое. Соответствующая разностная схема имеет вид:


\begin{gathered}\frac{\widehat{w}_{i,j}^{n+1}-\widehat{w}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2}}{\tau\!\!\not{\phantom{|}}\,2}= \frac{\widehat{\nu}_{i+1,j}^{n+1}-2\, \widehat{w}_{i,j}^{n+1}+ \widehat{w}_{i-1,j}^{n+1}}{h_y^2}\,,\quad i,j=0,\pm1,\pm2,\ldots,\quad n=0,1,\ldots,N-1;\\[2pt] \widehat{w}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}= \widehat{\nu}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\quad i,j=0,\pm1,\pm2,\ldots,\quad n=0,1,\ldots,N-1, \end{gathered}
(9.25)

или в традиционной записи:


\frac{\tau}{2h_y^2}\,\widehat{w}_{i,j-1}^{n+1}-\left(1+ \frac{\tau}{h_y^2}\right)\! \widehat{w}_{i,j}^{n+1}+ \frac{\tau}{2h_y^2}\, \widehat{w}_{i,j+1}^{n+1}=-\widehat{\nu}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,2},\quad i,j=0,\pm1,\pm2,\ldots,\quad n=0,1,\ldots,N-1.

Для каждого значения / требуется решить трехдиагональную систему линейных алгебраических уравнений, так как каждое уравнение содержит три неизвестных значения \widehat{w}_{i,j-1}^{n+1},\, \widehat{w}_{i,j}^{n+1},\, \widehat{w}_{i,j+1}^{n+1}, а остальные значения берутся с (n+1/2)-го слоя. Иными словами, при таком переходе схема является неявной. Искомые значения на целом слое вычисляются методом прогонки по направлению y. Далее полагают \widehat{u}_{i,j}^{n+1}= \widehat{w}_{i,j}^{n+1},~ i,j=0,\pm1,\pm2,\ldots.


Можно показать, что схема (9.24)–(9.25) метода дробных шагов имеет первый порядок аппроксимации по \tau и второй по h и является абсолютно устойчивой. Диагональные элементы в системах (9.24),(9.25) преобладают, следовательно, процесс прогонки устойчив. К достоинствам метода, достигнутым благодаря введению промежуточного слоя, следует отнести расщепление исходной задачи на две более простые, решаемые с помощью алгоритма скалярной прогонки. К недостаткам метода относится сравнительно низкая точность из-за того, что на каждом дробном шаге аппроксимируется только один пространственный оператор.


Методика расчетов


1. Задать сетку D_h= \bigl\{x_{i},y_{j},t^n\bigr\}, где x_{i}=ih,~ y_{j}= jh,~ t^n=n\tau, i,j=0,\pm1,\pm2,\ldots, n=0,1,\ldots,N; T= N\tau. Величины шагов в силу абсолютной устойчивости метода выбрать из условия удовлетворения точности расчетов.


2. Вычислить \widehat{u}_{i,j}^{0}= \widehat{\nu}_{i,j}^{0}= \psi_{i,j},~ i,j=0,\pm1,\pm2,\ldots,, где \psi_{i,j}=\psi(x_{i},y_{j}). Положить n=0.


3. Найти решение на (n+1)-м слое.


Первый этап. Перейти на дополнительный промежуточный (n+1/2)-й слой с шагом \tau\!\!\not{\phantom{|}}\,2\colon для каждого j решить трехдиагональную систему линейных алгебраических уравнений с учетом п.2:


\frac{\tau}{2h_x^2}\,\widehat{\nu}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2}-\left(1+ \frac{\tau}{h_x^2}\right)\!\widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2}+ \frac{\tau}{2 h_x^2}\, \widehat{\nu}_{i+1,j}^{n+1\!\not{\phantom{|}}\,\,2}=-\widehat{u}_{i,j}^{n},\quad i,j=0,\pm1, \pm2, \ldots

Искомые значения \widehat{\nu}_{i-1,j}^{n+1\!\not{\phantom{|}}\,\,2},\, \widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},\, \widehat{\nu}_{i+1, j}^{n+ 1\!\not{\phantom{|}}\,\,2} на промежуточном слое вычисляются методом прогонки по направлению x.


Второй этап. Перейти на (n+1)-й временной слой с промежуточного (n+1/2)-го с шагом \tau\!\!\not{\phantom{|}}\,2\colon для каждого значения / решить трехдиагональную систему линейных алгебраических уравнений


\frac{\tau}{2h_y^2}\,\widehat{w}_{i,j-1}^{n+1}-\left(1+ \frac{\tau}{h_y^2}\right)\! \widehat{w}_{i,j}^{n+1}+ \frac{\tau}{2 h_y^2}\, \widehat{w}_{i,j+1}^{n+1}=-\widehat{\nu}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},\quad i,j=0,\pm1, \pm2, \ldots

Искомые значения \widehat{w}_{i,j-1}^{n+1},\, \widehat{w}_{i,j}^{n+1},\, \widehat{w}_{i,j+ 1}^{n+1} на целом слое вычисляются методом прогонки по направлению y.


Положить \widehat{u}_{i,j}^{n+1}= \widehat{w}_{i,j}^{n+1},~ i,j=0,\pm1,\pm2,\ldots.


4. Если n=N, вычисления завершить. Иначе положить n=n+1 и перейти к пункту 3.


Замечания


1. В приведенной методике для простоты можно сразу положить \widehat{\nu}_{i, j}^{n+ 1\!\not{\phantom{|}}\,\,2}= \widehat{u}_{i,j}^{n+1\!\not{\phantom{|}}\,\,2},~ \widehat{w}_{i, j}^{n+1}= \widehat{u}_{i,j}^{n+1}.


2. При применении метода дробных шагов к решению уравнения с тремя пространственными переменными u_{t}= u_{xx}+ u_{yy}+ u_{zz} используются два промежуточных слоя: t=t^{n+1\!\not{\phantom{|}}\,\,3} и t=t^{n+2\!\not{\phantom{|}}\,\,3}. Тогда на дробных шагах реализуется следующая аппроксимация дифференциальных операторов:


\begin{aligned}&\frac{\widehat{\nu}_{i,j,k}^{n+ 1\!\not{\phantom{|}}\,\,3}-\widehat{\nu}_{i,j,k}^{n} }{\tau\!\!\not{\phantom{|}}\,3}= \frac{\widehat{\nu}_{i+1,j,k}^{n+ 1\!\not{\phantom{|}}\,\,3}-2\, \widehat{\nu}_{i,j,k}^{n+ 1\!\not{\phantom{|}}\,\,3}+ \widehat{\nu}_{i-1,j,k}^{n+ 1\!\not{\phantom{|}}\,\,3}}{h_x^2}\,, & \quad & \widehat{\nu}_{i,j,k}^{n}= \widehat{u}_{i,j,k}^{n},\\[4pt] &\frac{\widehat{w}_{i,j}^{n+ 2\!\not{\phantom{|}}\,\,3}-\widehat{w}_{i,j}^{n} }{\tau\!\!\not{\phantom{|}}\,3}= \frac{\widehat{\nu}_{i,j+1,k}^{n+ 2\!\not{\phantom{|}}\,\,3}-2\, \widehat{w}_{i,j}^{n+ 2\!\not{\phantom{|}}\,\,3}+ \widehat{\nu}_{i,j-1,k}^{n+ 2\!\not{\phantom{|}}\,\,3}}{h_y^2}\,, & \quad & \widehat{w}_{i,j}^{n+ 1\!\not{\phantom{|}}\,\,3}= \widehat{\nu}_{i,j,k}^{n+ 1\!\not{\phantom{|}}\,\,3},\\[4pt] &\frac{\widehat{u}_{i,j,k}^{n+1}-\widehat{u}_{i,j,k}^{n+ 2\!\not{\phantom{|}}\,\,3} }{\tau\!\!\not{\phantom{|}}\,3}= \frac{\widehat{u}_{i,j,k+1}^{n+1}-2\, \widehat{u}_{i,j,k}^{n+1}+ \widehat{\nu}_{i,j,k-1}^{n+1}}{h_z^2}\,, & \quad & \widehat{u}_{i,j,k}^{n+ 2\!\not{\phantom{|}}\,\,3}= \widehat{u}_{i,j,k}^{n+ 2\!\not{\phantom{|}}\,\,3}. \end{aligned}

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

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

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


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

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