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

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

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

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


Разностные схемы решения уравнений в частных производных 2-го порядка

Разностные схемы решения уравнений
в частных производных 2-го порядка


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


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


Пример 8.3. Построить явную разностную схему для решения задачи теплопроводности в стержне, начальная температура которого равна нулю. Пусть температура левого конца (x=0) фиксирована, а на правом конце (x=L) происходит теплообмен с окружающей средой, так что тепловой поток пропорционален разности температур конца стержня и среды (определяется функцией g(t)).


Решается третья начально-краевая задача:
u_t=u_{xx},~~ 0<x<L,~ 0<t<T,
u(x,0)=0,~~ 0 \leqslant x \leqslant L, (начальное условие)
u(0,t)=1,~~ 0<t<T, (краевое условие на левой границе)
u_x(L,t)=-\bigl[u(L,t)-g(t)\bigr],~~ 0<t \leqslant T, (краевое условие на правой границе).

Сравнивая с общей постановкой задачи, получаем (L,\,T — заданные числа)


a^2=1,\quad \psi(x)\equiv 0,,\quad \alpha_0=0,\quad \beta=1,\quad \varphi_1(x)=1,\quad \alpha_1=1,\quad \beta_1=1,\quad \varphi_2(t)=g(t).

Для аппроксимации производных в решаемом уравнении будем использовать (8.27),(8.26): \frac{\widehat{u}_{i}^{n+1}-\widehat{u}_{i}^{n}}{\tau}= \frac{\widehat{u}_{i+1}^{n}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i-1}^{n}}{h^2}. Отсюда


\widehat{u}_{i}^{n+1}= \frac{\tau}{h^2}\,\widehat{u}_{i-1}^{n}+ \left(1-\frac{2\tau}{h^2}\right)\!\widehat{u}_{i}^{n}+ \frac{\tau}{h^2}\,\widehat{u}_{i+1}^{n},\quad i=1,2,\ldots,I-1.
(8.41)

Эта формула при фиксированном значении i выражает решение на (n+1)-м временном слое через решение на n-м слое. Ей соответствует четырехточечный шаблон, изображенный на рис. 8.5, в.


При записи в разностной форме начальное условие u(x,0)=0,~ 0 \leqslant x \leqslant L, запишется в следующем виде:


\widehat{u}_{i}^{0}=0,\quad i=0,1,2,\ldots,I
(8.42)

краевые условия на левом конце (x=0)\colon


\widehat{u}_{0}^{n}=1,\quad n=1,2,\ldots,N,
(8.43)

краевые условия на правом конце (x=L) определяются левой разностью по формуле (8.21): \frac{\widehat{u}_{I}^{n}-\widehat{u}_{I-1}^{n}}{h}=-\bigl[\widehat{u}_{I}^{n}-g(n\,\tau)\bigr]. Отсюда


\widehat{u}_{I}^{n}= \frac{\widehat{u}_{I-1}^{n}-h\cdot g(n\cdot \tau)}{1+h},\quad n=1,2,\ldots,N.
(8.44)

Заметим, что соотношения (8.42),(8.43) точно аппроксимируют начальное и соответствующее краевое условие. Соотношения (8.41)–(8.44) образуют явную двухслойную разностную схему, имеющую первый порядок аппроксимации по \tau и второй — по h.


Замечания


1. Данная схема условно устойчива (устойчива при выполнении условия \frac{\tau}{h^2} \leqslant \frac{1}{2}, т.е. при \tau \leqslant \frac{h^2}{2}). Тем самым накладываются серьезные ограничения на выбор шага по времени. Например, если шаг сетки по переменной x выбрать равным h=0,\!1, то шаг сетки по времени не должен быть больше, чем \tau=0,\!5\cdot 0,\!1^2= 0,\!005 (это означает, что отрезок [0;1] будет пройден за 200 шагов). Наилучшим с точки зрения минимизации ошибок является соотношение \tau= \frac{h^2}{6}.


2. Разностная схема для уравнения u_t= a^2(x,t)u_{xx} строится аналогично и является устойчивой при \frac{\tau}{h^2} \leqslant \frac{1}{2A}, где A= \max_{(x,t)\in D} a^2(x,t).


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


1. Задать значения шагов сетки A,\,\tau так, что l\,h=L,~ N\,\tau=T, где I,\,N — целые положительные числа, причем значение шага по времени выбрать из условия устойчивости \tau \leqslant \frac{h^2}{2}. Вычислить \widehat{u}_{i }^{0}=0,~ i=0,1,\ldots,I; \widehat{u}_{0}^{n}=1,~ n=1,\ldots,N. Положить n=0.


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


\widehat{u}_{i}^{n}= \frac{\tau}{h^2}\,\widehat{u}_{i-1}^{n}+ \left(1-\frac{2\tau}{h^2}\right)\!\widehat{u}_{i}^{n}+ \frac{\tau}{h^2}\,\widehat{u}_{i+1}^{n},\quad i=1,\ldots,I-1;\quad \widehat{u}_{I}^{n}= \frac{\widehat{u}_{I-1}^{n}-h\,g(n\,\tau)}{1+h}\,.

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


Разностная схема на сетке 1

Пример 8.4. Пользуясь явной схемой, используемой в примере 8.3, найти приближенное решение первой начально-краевой задачи вида:

u_t=u_{xx},~ 0<x<1,~0<t<0,\!01,

u(x,0)= 4x(1-x),~ 0 \leqslant x \leqslant 1, (начальное условие)

u(0,t)=0,~ 0<t \leqslant 0,\!01, (краевое условие на левой границе)

u(1,t)=0,~ 0<t \leqslant 0,\!01, (краевое условие на правой границе).


Сравнивая с общей постановкой задачи (см. разд. 8.2), получаем


a^2=1,\quad \psi(x)=4x(1-x),\quad \varphi_1(t)\equiv 0,\quad \varphi_2(t)\equiv 0,\quad L=1,\quad T=0,\!01.

1. Зададим шаг по пространству: h=0,\!01, следовательно, I=\frac{L}{0,\!1}= \frac{1}{0,\!1}=10. Выберем шаг по времени из условия обеспечения устойчивости и минимизации ошибок \tau=\frac{h^2}{6}=\frac{0,\!01}{6}=\frac{1}{600}, поэтому N= \frac{T}{\tau}= \frac{0,\!01}{1\!\!\not{\phantom{|}}\,600}=6. Таким образом, получаем


x_{i}=0,\!1i,\quad i=0,1,2,\ldots,10;\qquad t^n= \frac{n}{600},\quad n=0,1,2,\ldots,6.

Вычислим значения решения на нулевом слое:


\widehat{u}_{i}^{0}= 4x_{i}(1-x_{i})= 4\cdot 0,\!1\cdot i(1-0,\!1i),\quad i=0,1,2,\ldots,10~ (I=10).

Краевые условия на левой и правой границах записываются в виде:


\widehat{u}_{0}^{n}=0,\quad n=1,2,\ldots,N=6;\qquad \widehat{u}_{I}^{n}=0,\quad n=1,2,\ldots,6~ (N=6).

2,3. Поскольку \frac{\tau}{h^2}= \frac{1\!\!\not{\phantom{|}}\,600}{0,\!01}= \frac{1}{6} решение на следующих слоях при n=0,1,\ldots,5 находится по формуле, следующей из (8.41):


\widehat{u}_{i}^{n+1}= \frac{1}{6}\cdot \widehat{u}_{i-1}^{n}+ \frac{2}{3}\cdot \widehat{u}_{i}^{n}+ \frac{1}{6}\cdot \widehat{u}_{i+1}^{n},\quad i=1,2,\ldots,9;

и занесено в табл. 8.1.

\begin{aligned}\mathit{Table~8.1}~&\\ \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|c|} \hline n \setminus i& 0& 1& 2& 3& 4& 5& 6& 7& 8& 9& 10 \\\hline 0& 0& 0,\!360& 0,\!640& 0,\!840& 0,\!960& 1,\!000& 0,\!960& 0,\!840& 0,\!640& 0,\!360& 0 \\\hline <br />1& 0& 0,\!347& 0,\!627& 0,\!827& 0,\!947& 0,\!987& 0,\!947& 0,\!827& 0,\!627& 0,\!347& 0 \\\hline 2& 0& 0,\!336& 0,\!613& 0,\!813& 0,\!933& 0,\!973& 0,\!933& 0,\!813& 0,\!613& 0,\!336& 0 \\\hline 3& 0& 0,\!326& 0,\!600& 0,\!800& 0,\!920& 0,\!960& 0,\!920& 0,\!800& 0,\!600& 0,\!326& 0 \\\hline 4& 0& 0,\!317& 0,\!588& 0,\!787& 0,\!907& 0,\!947& 0,\!907& 0,\!787& 0,\!588& 0,\!317& 0 \\\hline 5& 0& 0,\!309& 0,\!576& 0,\!774& 0,\!894& 0,\!934& 0,\!894& 0,\!774& 0,\!576& 0,\!309& 0 \\\hline 6& 0& 0,\!302& 0,\!564& 0,\!761& 0,\!881& 0,\!921& 0,\!881& 0,\!761& 0,\!564& 0,\!302& 0 \\\hline \end{array}& \end{aligned}

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


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


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


Пример 8.5. Построить неявную разностную схему решения задачи теплопроводности в стержне, начальная температура которого равна 1, а температура левого и правого концов равна нулю.


Данным физическим условиям соответствует первая начально-краевая задача:


u_t=u_{xx},~~ 0&lt;x&lt;L,~ 0&lt;t&lt;T,
u(x,0)=1,~~ 0 \leqslant x \leqslant L, (начальное условие)
u(0,t)=0,~~ 0&lt;t \leqslant T, (краевое условие на левой границе)
u(L,t)=0,~~ 0&lt;t \leqslant T (краевое условие на правой границе).

Сравнивая с общей постановкой задачи, получаем a^2=1,~ \psi(x)\equiv 1~ \varphi_1(t)= \varphi_2(t)\equiv 0; L,\,T— заданные числа.


Для аппроксимации производных в решаемом уравнении будем использовать формулу (8.27) и формулу (8.26), записанную на (n+1)-м временном слое:


\frac{\widehat{u}_{i}^{n+1}-\widehat{u}_{i}^{n}}{\tau}= \frac{\widehat{u}_{i+1}^{n+1}-2\widehat{u}_{i}^{n+1}+ \widehat{u}_{i-1}^{n+1}}{h^2}\,.

Отсюда

-\frac{\tau}{h^2}\cdot \widehat{u}_{i-1}^{n+1}+ \left(1+\frac{2\tau}{h^2}\right)\cdot \widehat{u}_{i}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{i+1}^{n+1}= \widehat{u}_{i}^{n},\quad i=1,2,\ldots,I-1.
(8.45)

Эта формула выражает решение на (n+1)-м временном слое через решение на n-м слое. Ей соответствует четырехточечный шаблон, изображенный на рис. 8.8,а.


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

При записи в разностной форме начальное условие u(x,0)=1,~ 0 \leqslant x \leqslant L, запишется в следующем виде:


\widehat{u}_{i}^{0}=1,\quad i=0,1,2,\ldots,I,
(8.46)

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

\begin{aligned}&\widehat{u}_{0}^{n}=0, & & n=1,2,\ldots,N,\\ &\widehat{u}_{I}^{n}=0, & & n=1,2,\ldots,N. \end{aligned}
(8.47)

Соотношения (8.45)–(8.47) образуют неявную двухслойную разностную схему, имеющую первый порядок аппроксимации по \tau и второй — по h. Она является абсолютно устойчивой и сходящейся.


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


1. Задать значения шагов сетки h,\,\tau из условия обеспечения требуемой точности решения так, что I\,h=L,~ N\,\tau=T,, где I,\,N — целые положительные числа. Вычислить \widehat{u}_{i}^{0}=1,~ i=0,1,\ldots,I; \widehat{u}_{0}^{n}=0,~ n=1,\ldots,N; \widehat{u}_{I}^{n}=0,~ n=1,\ldots,N. Положить n=0.


2. Сформируем систему линейных алгебраических уравнений для некоторого временного слоя, записывая уравнение (8.45) при i=1,2,\ldots,I-1 (рис. 8.9):


\begin{aligned}&-\frac{\tau}{h^2}\cdot \widehat{u}_{0}^{n+1}+ \left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{1}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{2}^{n+1}= \widehat{u}_{1}^{n},\\[2pt] &-\frac{\tau}{h^2}\cdot \widehat{u}_{1}^{n+1}+ \left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{2}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{3}^{n+1}= \widehat{u}_{2}^{n},\\[-2pt] &\quad\vdots\\[-4pt] &-\frac{\tau}{h^2}\cdot \widehat{u}_{I-2}^{n+1}+ \left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{I-1}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{I}^{n+1}= \widehat{u}_{I-1}^{n}. \end{aligned}

Согласно п.1, \widehat{u}_{0}^{n}=1,~ n=1,\ldots,N;~ \widehat{u}_{1}^{n}=0,~ n=1,\ldots,N. Поэтому для нахождения решения на следующем временном слое требуется решить трехдиагональную систему (I-1) линейных алгебраических уравнений с (I-1) неизвестными \widehat{u}_{1}^{n+1},\ldots,\widehat{u}_{I-1}^{n+1}\colon


\begin{gathered}\left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{1}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{2}^{n+1}= \widehat{u}_{1}^{n},\\ -\frac{\tau}{h^2}\cdot \widehat{u}_{1}^{n+1}+ \left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{2}^{n+1}-\frac{\tau}{h^2}\cdot \widehat{u}_{3}^{n+1}= \widehat{u}_{2}^{n},\\[-2pt] \vdots\\[-4pt]-\frac{\tau}{h^2}\cdot \widehat{u}_{I-2}^{n+1}+ \left(1+ \frac{2\tau}{h^2}\right)\cdot \widehat{u}_{I-1}^{n+1}= \widehat{u}_{I-1}^{n}. \end{gathered}

Методы решения таких систем рассматривались ранее.


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


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


Разностная схема на сетке 2

Получим еще одну неявную схему. Для аппроксимации производной щ применим (8.27), а для аппроксимации и^ взвешенное среднее выражения (8.26), записанного для узлов (x_{i}, t^n),\, (x_{i},t^{n+1}) на n-м и (n+1)-м временных слоях (см. рис. 8.8,б):


\widehat{u}_{xx}= \frac{\widehat{u}_{i+1}^{n}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i-1}^{n}}{h^2}\,,\qquad \widehat{u}_{xx}= \frac{\widehat{u}_{i+1}^{n+1}-2\widehat{u}_{i}^{n+1}+ \widehat{u}_{i-1}^{n+1}}{h^2}\,,

то есть u_{xx}= \lambda\,\frac{\widehat{u}_{i+1}^{n+1}-2\widehat{u}_{i}^{n+1}+ \widehat{u}_{i-1}^{n+1}}{h^2}+ (1-\lambda) \frac{\widehat{u}_{i+1}^{n}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i-1}^{n}}{h^2}, где \lambda\in [0;1]. При \lambda= \frac{1}{2} получается обычное среднее, при \lambda= \frac{3}{4} одна из разностных производных берется с весом \frac{3}{4}, а другая — с весом \frac{1}{4}. В результате получаем разностную схему:


\begin{gathered}\frac{\widehat{u}_{i}^{n+1}-\widehat{u}_{i}^{n}}{\tau}= \lambda\,\frac{\widehat{u}_{i+ 1}^{n+1}-2\widehat{u}_{i}^{n+1}+ \widehat{u}_{i-1}^{n+1}}{h^2}+ (1-\lambda) \frac{\widehat{u}_{i+ 1}^{n}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i-1}^{n}}{h^2}\,,\quad i=1,2,\ldots,I-1,\\[2pt] \widehat{u}_{i}^{0}=1,~ i=0,1,\ldots,I;\quad \widehat{u}_{0}^{n}=1,~ n=1,2,\ldots,N;\quad \widehat{u}_{I}^{n}=0,\quad n=1,2,\ldots,N. \end{gathered}

Перенесем в первых (I-1) уравнениях все неизвестные в левую часть. Тогда получаем


\begin{gathered}-\lambda\cdot \frac{\tau}{h^2}\cdot \widehat{u}_{i-1}^{n+1}+ \left(1+2 \lambda \frac{\tau}{h^2}\right)\cdot \widehat{u}_{i}^{n+1}-\lambda\cdot \frac{\tau}{h^2}\cdot \widehat{u}_{i+1}^{n+1}=\\ =(1-\lambda)\cdot \frac{\tau}{h^2}\cdot \widehat{u}_{i-1}^{n}+ \left(1+2 (1-\lambda) \frac{\tau}{h^2}\right)\cdot \widehat{u}_{i}^{n}-(1-\lambda)\cdot \frac{\tau}{h^2}\cdot \widehat{u}_{i+1}^{n},\quad i=1,\ldots,I-1;\\[2pt] \widehat{u}_{i}^{0}=1,~ i=0,1,\ldots,I;\quad \widehat{u}_{0}^{n}=0,~ n=1,2,\ldots,N;\quad \widehat{u}_{I}^{n}=1,~ n=1,2,\ldots,N. \end{gathered}
(8.48)

Эта схема абсолютно устойчива и сходится. Для нахождения решения на следующем слое, как и ранее, требуется решать трехдиагональную систему (I-1) линейных алгебраических уравнений. При \lambda= \frac{1}{2} схема называется схемой Кранка-Николсон. При \lambda=0 из (8.48) следует явная схема (8.41), а при \lambda=1 — неявная схема (8.45). Шеститочечный шаблон, соответствующий схеме (8.48), изображен на рис. 8.8,б. При \lambda= \frac{1}{2} схема имеет второй порядок аппроксимации по \tau и по h, а при \lambda\ne \frac{1}{2} — первый порядок аппроксимации по \tau и второй — по h.


Методику вычислений, совпадающую с уже описанной для схемы (8.45), продемонстрируем на примере.


1. Зададим h= 0,\!2;~ \tau= 0,\!08;~ L= 1;~ T= 0,\!8. Тогда I\,0,\!2= 1,~ N\,0,\!08 = 0,\!8, то есть I=5,~ N=10. Тогда начальные и краевые условия примут вид:


\widehat{u}_{i}^{0}=1,~ i=0,1,\ldots,5;\quad \widehat{u}_{0}^{n}=0,~ n=1,2,\ldots,10;\quad \widehat{u}_{5}^{n}=0,~ n=1,2,\ldots,10.

Положим n=0,~ \lambda= \frac{1}{2}. Вычислим \frac{\tau}{h^2}= \frac{0,\!08}{0,\!04}=2.


2. Запишем систему (8.48) , принимая i=1,2,3,4~ (I-1=4) с учетом пункта 1:


\begin{aligned}&-\widehat{u}_{0}^{1}+3\widehat{u}_{1}^{1}-\widehat{u}_{2}^{1}= \widehat{u}_{0}^{0}-\widehat{u}_{1}^{0}+\widehat{u}_{2}^{0},\\[2pt] &-\widehat{u}_{1}^{1}+3\widehat{u}_{2}^{1}-\widehat{u}_{3}^{1}= \widehat{u}_{1}^{0}-\widehat{u}_{3}^{0}+\widehat{u}_{3}^{0},\\[2pt] &-\widehat{u}_{2}^{1}+3\widehat{u}_{3}^{1}-\widehat{u}_{4}^{1}= \widehat{u}_{2}^{0}-\widehat{u}_{3}^{0}+\widehat{u}_{4}^{0},\\[2pt] &-\widehat{u}_{3}^{1}+3\widehat{u}_{4}^{1}-\widehat{u}_{5}^{1}= \widehat{u}_{3}^{0}-\widehat{u}_{4}^{0}+\widehat{u}_{5}^{0};\\[4pt] &\widehat{u}_{i}^{0}=1,~ i=0,1,\ldots,5;\quad \widehat{u}_{0}^{1}=0,\quad \widehat{u}_{5}^{1}=0. \end{aligned}

Тогда получаем \begin{cases}3\,\widehat{u}_{1}^{1}-\widehat{u}_{2}^{1}=1,\\-\widehat{u}_{1}^{1}+ 3\,\widehat{u}_{2}^{1}-\widehat{u}_{3}^{1}=1,\\-\widehat{u}_{2}^{1}+ 3\,\widehat{u}_{3}^{1}-\widehat{u}_{4}^{1}=1,\\-\widehat{u}_{3}^{1}+ 3\,\widehat{u}_{4}^{1}=1; \end{cases} или в матричной форме \begin{pmatrix}3&-1&0&0\\-1&3&-1&0\\ 0&-1&3&-1\\ 0&0&-1&3 \end{pmatrix}\!\cdot\! \begin{pmatrix}\widehat{u}_{1}^{1}\\ \widehat{u}_{2}^{1}\\ \widehat{u}_{3}^{1}\\ \widehat{u}_{4}^{1} \end{pmatrix}= \begin{pmatrix} 1\\1\\1\\1 \end{pmatrix}\!. В результате имеем \begin{cases}\widehat{u}_{1}^{1}=0,\!6;\\ \widehat{u}_{2}^{1}=0,\!8;\\ \widehat{u}_{3}^{1}=0,\!8;\\ \widehat{u}_{4}^{1}=0,\!6. \end{cases}


3. Далее, полагая n=1, можно сделать следующий шаг по времени, решая новую систему алгебраических уравнений. Аналогичный процесс продолжается до достижения значения N=10.




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


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


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


Пример 8.6. Построить явную и неявную разностные схемы для решения задачи, в которой имеется струна длиной L, натянутая между двумя точками оси Ox, точкой x=0 и точкой x=L. Концы струны закреплены, начальное смещение струны описывается функцией f(x), а начальная скорость — функцией g(x).


Решается первая начально-краевая задача:
u_{tt}=u_{xx},~~ 0&lt;x&lt;T,~ 0&lt;t&lt;T,
u(x,0)=f(x),~~ 0&lt;x&lt;L, (начальное условие)
u_t(x,0)=g(x),~~ 0&lt;x&lt;L, (начальное условие)
u(0,t)=0,~~ 0 \leqslant t \leqslant T, (краевое условие на левой границе)
u(L,t)=0,~~ 0 \leqslant t \leqslant T, (краевое условие на правой границе).

Сравнивая с общей постановкой задачи, получаем c^2=1,~ \psi_1(x)=f(x), \psi_2(x)= g(x), \varphi_1(t)= \varphi_2(t)\equiv 0; L,\,T — заданные числа.


Для аппроксимации производных в решаемом уравнении будем использовать (8.30) и (8.26):


\frac{\widehat{u}_{i}^{n+1}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i}^{n-1}}{\tau^2}= \frac{\widehat{u}_{i+1}^{n}-2\widehat{u}_{i}^{n}+\widehat{u}_{i-1}^{n}}{h^2}\,.

Отсюда


\widehat{u}_{i}^{n+1}= \frac{\tau^2}{h^2}\cdot \widehat{u}_{i-1}^{n}+ 2\! \left(1-\frac{\tau^2}{h^2}\right)\cdot \widehat{u}_{i}^{n}+ \frac{\tau^2}{h^2}\cdot \widehat{u}_{i+1}^{n}-\widehat{u}_{i}^{n-1},\quad i=1,2,\ldots,I-1.
(8.49)

Эта формула выражает решение на (n+1)-м временном слое через решение на n-м и (n-1)-м слоях. В отличие от схем для уравнения теплопроводности, в которых использовались только два временных слоя (n-й и (n+1)-й) , здесь требуется использовать три слоя: (n-1)-й, n-й и (n+1)-й. Ей соответствует пятиточечный шаблон, изображенный на рис. 8.10,а , называемый "крест".


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

При записи в разностной форме начальное условие u(x,0)= f(x),~ 0&lt;x&lt;L, запишется в следующем виде:


\widehat{u}_{i}^{0}= f(i\,h),\quad i=1,2,\ldots,I-1,
(8.50)

начальное условие u_t(x,0)=g(x),~ 0&lt;x&lt;L, с применением (8.27) при n=0 в виде \frac{\widehat{u}_{i}^{1}-\widehat{u}_{i}^{0}}{\tau}= g(i\,h), i=1,2,\ldots,I-1, или, используя предыдущее соотношение,


\widehat{u}_{i}^{1}= f(i\,h)+\tau\cdot g(i\,h),\quad i=1,2,\ldots,I-1,
(8.51)

краевые условия представляются в форме:


\widehat{u}_{0}^{n}=0,\quad \widehat{u}_{I}^{n}=0,\quad n=0,1,\ldots,N.
(8.52)

Соотношения (8.49)–(8.52) образуют явную трехслойную разностную схему. Эта схема имеет первый порядок аппроксимации по t и второй — по x, так как соотношение (8.51) аппроксимирует дифференциальное начальное условие с первым порядком.


Замечания


1. Данная схема условно устойчива (устойчива при выполнении условия \frac{\tau}{h}\leqslant 1.


2. Разностная схема для уравнения u_{tt}= c^2(x,t)u_{xx} строится аналогично и является устойчивой при \frac{\tau^2}{h^2}&lt; \frac{1}{C}, где C= \max_{(x,t)\in D} c^2(x,t).


3. Простейшая замена (8.51) имеет первый порядок аппроксимации. Поскольку разностная схема (8.49) аппроксимирует дифференциальное уравнение со вторым порядком, желательно, чтобы разностное начальное условие также имело второй порядок аппроксимации. Для этого запишем аналог формулы Тейлора (8.17) в направлении t при k=2\colon


u(x,t^0+\tau)= u(x,t^0)+ u_t(x,t^0)\tau+ u_{tt}(x,t^0)\frac{\tau^2}{2}+ u_{ttt}(x,\xi)\frac{\tau^3}{6}\,,\quad \xi\in (t^0,t^0+\tau),~ 0&lt;x&lt;L.

Если функция f(x) имеет ограниченную вторую производную, то в силу исходного дифференциального уравнения u_{tt}= u_{xx} получаем u_{tt}(x,t^0)= u_{xx}(x,t^0)= f''(x). Из предыдущего соотношения получаем


u_{t}(x,t^0)= \frac{u(x,t^0+\tau)-u(x,t^0)}{\tau}+ f''(x)\frac{\tau^2}{2}-u_{ttt}(x,\xi)\frac{\tau^3}{6}

или

u_{t}(x,t^0)\cong \frac{u(x,t^0+\tau)-u(x,t^0)}{\tau}+ f''(x)\frac{\tau^2}{2}\,,\quad \left(\frac{\tau^2}{6}\,M_3\right)\!,

где M_3= \max_{t\in [t^0,t^0+\tau]} \bigl|u_{ttt}(x,t)\bigr|. В результате соответствующее начальное условие в разностной форме имеет вид \frac{\widehat{u}_{i}^{1}-\widehat{u}_{i}^{0}}{\tau}-f''(i\,h)\frac{\tau}{2}= g(i\,h) i=1,\ldots,I-1.


Так как \widehat{u}_{i}^{0}= f(i\,h),~ i=1,\ldots,I-1, то


\widehat{u}_{i}^{1}= f(ih)+ \tau\cdot g(ih)+ f''(ih)\cdot \frac{\tau^2}{2}\,,\quad i=1,2,\ldots,I-1.
(8.53)

Если в разностной схеме используется (8.51), схема имеет первый порядок аппроксимации по \tau и второй — по h, а если (8.53), то имеет второй порядок аппроксимации по \tau и по h.


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


1. Задать значения шагов сетки h,\,\tau так, что I\,h=L,~ N\,\tau=T, где I,\,N — целые положительные числа, причем значение шага по времени выбрать из условия устойчивости \tau \leqslant h. Вычислить


\begin{gathered}\widehat{u}_{i}^{0}= f(ih),\quad i=1,2,\ldots,I-1;\qquad \widehat{u}_{i}^{1}= f(ih)+\tau g(ih),\quad i=1,2,\ldots,I-1;\\[2pt] \widehat{u}_{i}^{0}=0,\quad \widehat{u}_{I}^{n}=0,\quad n=0,1,\ldots,N. \end{gathered}

Приведенные выражения дают решение задачи для первых двух слоев сетки (рис. 8.8). Положить n=1.


Разностная схема на сетке 3

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


\widehat{u}_{i}^{n+1}= \frac{\tau^2}{h^2}\cdot \widehat{u}_{i-1}^{n}+ 2\! \left(1-\frac{\tau^2}{h^2}\right)\cdot \widehat{u}_{i}^{n}+ \frac{\tau^2}{h^2}\cdot \widehat{u}_{i+1}^{n}-\widehat{u}_{i}^{n-1},\quad i=1,2,\ldots,I-1.

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


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


Для иллюстрации общего подхода получим неявную схему для волнового уравнения. При аппроксимации производных в уравнении будем использовать (8.30) и формулу (8.26), записанную на (n-1)-м, n-м и (n+1)-м слоях:


\begin{aligned}\frac{\widehat{u}_{i}^{n+1}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i}^{n-1}}{\tau^2}&= \lambda\cdot \frac{\widehat{u}_{i+1}^{n+1}-2\widehat{u}_{i}^{n+ 1}+\widehat{u}_{i-1}^{n+1}}{h^2}+ (1-2 \lambda)\cdot \frac{\widehat{u}_{i+ 1}^{n}-2\widehat{u}_{i}^{n}+ \widehat{u}_{i-1}^{n}}{h^2}\,+\\ &\quad +\,\lambda\cdot \frac{\widehat{u}_{i+1}^{n-1}-2\widehat{u}_{i}^{n-1}+\widehat{u}_{i-1}^{n-1}}{\tau^2}\,\quad i=1,2,\ldots,I-1,\end{aligned}
(8.54)

где \lambda — параметр, 0 \leqslant \lambda \leqslant \frac{1}{2}. Значения приближенного решения на нулевом и первом слоях вычисляются по формулам (8.50)–(8.52), как и в явной схеме. На остальных слоях схема представляет собой линейную систему алгебраических уравнений с трехдиагональной матрицей, в которой диагональные элементы преобладают. Решение этой системы существует, единственно и вычисляется методом прогонки.


Формула (8.54) выражает три неизвестных значения решения на (n+1)-м временном слое через решения на n-м и (n-1)-м слоях. В отличие от схем для уравнения теплопроводности, в которых использовались только два временных слоя (n-й и (n+1)-й) , здесь требуется использовать три слоя: (n-1)-й, n-й и (n+1)-й. Поэтому такая схема называется неявной трехслойной. Ей соответствует девятиточечный шаблон, изображенный на рис. 8.10,б. При \lambda=0 схема (8.50)-(8.52),(8.54) переходит в явную схему "крест".


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


Пример 8.7. Пользуясь явной схемой, найти приближенное решение первой начально-краевой задачи:
u_{tt}=u_{xx},~~ 0&lt;x&lt;\pi,~ 0&lt;t&lt;\frac{5\pi}{18}\,,
u(x,0)= x(\pi-x),~~ 0&lt;x&lt;\pi, (начальное условие)
u_t(x,0)=0,~~ 0&lt;x&lt;\pi (начальное условие)
u(0,t)=0,~~ 0 \leqslant t \leqslant \frac{5\pi}{18} (краевое условие на левой границе)
u(\pi,t)=0,~~ 0 \leqslant t \leqslant \frac{5\pi}{18} (краевое условие на правой границе).

Сравнивая с общей постановкой задачи, получаем


c^2=1,\quad \psi_1(x)=x(\pi-x),\quad \psi_2(x)\equiv0,\quad \varphi_1(t)= \varphi_2(t)\equiv 0;\quad L=\pi,\quad T=\frac{5\pi}{18}\,.

1. Зададим шаг по пространству, равным h=\frac{\pi}{18}, шаг по времени выберем из условия устойчивости: \tau \leqslant h. Положим \tau=h= \frac{\pi}{18}. Тогда


x_{i}=\frac{\pi\,i}{18}\,,\quad t^n= \frac{\pi\,n}{18}\,,\quad i=0,1,\ldots,I=\frac{L}{h}=18;\quad n=0,1,\ldots,N=\frac{T}{\tau}=5.

По формулам (8.50), (8.52), (8.53), учитывая, что f(x)=x(\pi-x), f''(x)=-2,~ g(x)\equiv 0, получаем значения приближенного решения на нулевом и первом слоях:


\begin{gathered}\widehat{u}_{i}^{0}= f(ih)= x_{i}(\pi-x_{i})= \frac{\pi\,i}{18}\! \left(\pi-\frac{\pi\,i}{18}\right)\!,\quad i=1,2,\ldots,17,\\[2pt] \widehat{u}_{0}^{n}=0,\quad \widehat{u}_{I}^{n}= 0,\quad n=0,\ldots,5,\\[2pt] \widehat{u}_{i}^{1}= f(ih)+ \tau\,g(ih)+ f''(ih)\frac{\tau^2}{2}= x_{i}(\pi-x_{i})+0-2\, \frac{\tau^2}{2}= \frac{\pi\,i}{18}\! \left(\pi-\frac{\pi\,i}{18}\right)-\left(\frac{\pi}{18} \right)^2\!,\quad i=1,2,\ldots,17. \end{gathered}

2,3. Поскольку \frac{\tau^2}{h^2}=1, решение на последующих слоях вычисляется по И формуле, следующей из (8.49):


\widehat{u}_{i}^{n+1}= \widehat{u}_{i-1}^{n}+ \widehat{u}_{i+1}^{n}-\widehat{u}_{i}^{n-1},\quad i=1,2,\ldots,17.

В табл. 8.2 приведены результаты расчетов при i=0,1,\ldots,9, так как решение u(x,t) симметрично относительно x=\frac{\pi}{2}~(i=9) при фиксированном t.


\begin{aligned}\mathit{Table~8.2}~&\\ \begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}\hline n \setminus i& 0& 1& 2& 3& 4& 5& 6& 7& 8& 9 \\\hline 0& 0& 0,\!518& 0,\!975& 1,\!371& 1,\!706& 1,\!980& 2,\!193& 2,\!346& 2,\!437& 2,\!467 \\\hline 1& 0& 0,\!487& 0,\!944& 1,\!340& 1,\!675& 1,\!950& 2,\!163& 2,\!315& 2,\!406& 2,\!437 \\\hline 2& 0& 0,\!426& 0,\!853& 1,\!249& 1,\!584& 1,\!858& 2,\!071& 2,\!224& 2,\!315& 2,\!346 \\\hline 3& 0& 0,\!366& 0,\!731& 1,\!097& 1,\!432& 1,\!706& 1,\!919& 2,\!071& 2,\!163& 2,\!193 \\\hline 4& 0& 0,\!305& 0,\!609& 0,\!914& 1,\!218& 1,\!493& 1,\!706& 1,\!858& 1,\!950& 1,\!980 \\\hline 5& 0& 0,\!244& 0,\!487& 0,\!731& 0,\!975& 1,\!218& 1,\!432& 1,\!584& 1,\!675& 1,\!706 \\\hline \end{array}& \end{aligned}

Замечание. В практике применения численных методов существует полезный прием, позволяющий на основе результатов расчетов судить о том, с какой точностью они получены. С этой целью используется правило Рунге, изложенное выше применительно к уточнению результатов численного дифференцирования и интегрирования. Аналогичный подход можно использовать для оценки точности решения задач математической физики. Опишем идею его реализации:

а) находится численное решение с шагом h по пространственной переменной x и шагом \tau по времени t;

б) расчет повторяется с шагами \frac{h}{2} и \frac{\tau}{2} с тем, чтобы образовывались слои, некоторые из которых совпадали бы со слоями, полученными ранее;

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




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


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


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


Пример 8.8. Построить разностную схему для решения задачи Дирихле (где \varphi(x,y) — заданная функция, L,\,M — заданные числа):

\begin{aligned}&u_{xx}+u_{yy}=0,~~ 0&lt;x&lt;L,~0&lt;y&lt;M,\\ &u(0,y)= \varphi(0,y),~~ 0&lt;y&lt;M,\\ &u(L,y)= \varphi(L,y),~~ 0&lt;y&lt;M,\\ &u(x,M)= \varphi(x,M),~~ 0 \leqslant x \leqslant L,\\ &u(x,0)= \varphi(x,0),~~ 0 \leqslant x \leqslant L, \end{aligned}


В данном случае задача решается в области П прямоугольной формы, краевые условия заданы на сторонах прямоугольника. В плоскости Oxy построим сетку D_h= \bigl\{(x_{i}, y_{i})\bigr\}, где x_{i}=ih_x,~ y_{j}=jh_{y}, i=0,1,\ldots,I; j=0,1,\ldots,J; h_x,\,h_y — величины шагов по x и y;~ I\,h_x=L, J\,h_y=M,~ I и J — целые положительные числа.


Заменим дифференциальные операторы в уравнении Лапласа по формулам, следующим из (8.25) при замене переменной t^n на переменную y_{j}\colon


\begin{aligned}& \widehat{u}_{xx}(x_{i},y_{j})= \frac{u(x_{i}+h_x,y_{j})-2u(x_{i},y_{j})+ u(x_{i}-h_x, y_{j})}{h_x^2}\,,\\[2pt] & \widehat{u}_{yy}(x_{i},y_{j})= \frac{u(x_{i},y_{j}+h_y)-2u(x_{i},y_{j})+ u(x_{i}, y_{j}-h_y)}{h_y^2}\,. \end{aligned}

Применяя краткую запись u(x_{i},y_{j})= u_{i,j}, получаем


\frac{\widehat{u}_{i+1,j}-2\widehat{u}_{i,j}+ \widehat{u}_{i-1,j}}{h_x^2}+ \frac{\widehat{u}_{i,j+1}-2\widehat{u}_{i,j}+\widehat{u}_{i,j-1}}{h_y^2}=0,\quad i=1,\ldots,I-1;~ j=1,\ldots,J-1,

или

\frac{h_y^2}{h_x^2}\,\widehat{u}_{i+1,j}+ \frac{h_y^2}{h_x^2}\,\widehat{u}_{i-1,j}+ \widehat{u}_{i,j+1}+ \widehat{u}_{i,j-1}-2\! \left(1+ \frac{h_y^2}{h_x^2}\right)\! \widehat{u}_{i,j}=0,\quad i=1,\ldots,I-1;~ j=1,\ldots,J-1.
(8.55)

Полагая \varphi_{i,j}= \varphi(ih_x,jh_y), краевые условия можно записать в виде:


\begin{aligned}& \widehat{u}_{0,j}= \varphi_{0,j},\quad j=1,\ldots,J-1;\\ & \widehat{u}_{I,j}= \varphi_{I,j},\quad j=1,\ldots,J-1;\\ & \widehat{u}_{i,J}= \varphi_{i,J},\quad i=0,1,\ldots,I;\\ & \widehat{u}_{i,0}= \varphi_{i,0},\quad i=0,1,\ldots,I. \end{aligned}
(8.56)

Для нахождения искомого решения требуется решить систему (I-1)\times (J-1) линейных алгебраических уравнений. Схема (8.55),(8.56) имеет второй порядок аппроксимации по \tau и по h, является абсолютно устойчивой и сходится. Ей соответствует пятиточечный шаблон, изображенный на рис. 8.12. При h_x= h_y, т.е. одинаковых шагах в горизонтальном и вертикальном направлениях, получаем


\widehat{u}_{i+1,j}+ \widehat{u}_{i-1,j}+ \widehat{u}_{i,j+1}+ \widehat{u}_{i,j-1}-4\,\widehat{u}_{i,j}=0,\quad i=1,\ldots,I-1;~ j=1,\ldots,J-1.
(8.57)

Разрешив относительно \widehat{u}_{i,j}, имеем


\widehat{u}_{i,j}= \frac{1}{4}\bigl(\widehat{u}_{i+1,j}+ \widehat{u}_{i-1,j}+ \widehat{u}_{i,j+1}+ \widehat{u}_{i,j-1}\bigr),\quad i=1,\ldots,I-1;~ j=1,\ldots,J-1.
(8.58)

Это означает, что решение \widehat{u}_{i,j} аппроксимируется средним значением по четырем соседним узлам. Как правило, для решения системы (8.58) применяются различные итерационные методы, например, метод простой итерации или метод Зейделя.


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

Рассмотрим частный случай поставленной задачи при L=1,~ M=1,~ I=J=3, \varphi(x,y)= \begin{cases}\sin\pi x,& y=0,~ 0 \leqslant x \leqslant 1,\\ 0,& \text{na ostalnih granicah}~\Omega. \end{cases} Так как h_x= h_y= \frac{1}{3}, то запишем систему (8.57) для четырех внутренних точек (i=1,2;~ j=1,2)\colon


\begin{aligned}& \widehat{u}_{2,1}+ \widehat{u}_{0,1}+ \widehat{u}_{1,2}+ \widehat{u}_{1,0}-4\,\widehat{u}_{1,1}=0,\quad (i=1,~ j=1)\\[2pt] & \widehat{u}_{2,2}+ \widehat{u}_{0,2}+ \widehat{u}_{1,3}+ \widehat{u}_{1,1}-4\,\widehat{u}_{1,2}=0,\quad (i=1,~ j=2)\\[2pt] & \widehat{u}_{3,1}+ \widehat{u}_{1,1}+ \widehat{u}_{2,2}+ \widehat{u}_{2,0}-4\,\widehat{u}_{2,1}=0,\quad (i=2,~ j=1)\\[2pt] & \widehat{u}_{3,2}+ \widehat{u}_{1,2}+ \widehat{u}_{2,3}+ \widehat{u}_{2,1}-4\,\widehat{u}_{2,2}=0.\quad (i=2,~ j=2) \end{aligned}

Краевые условия (8.56) запишутся в форме


\begin{aligned}& \widehat{u}_{0,j}=0,~~ j=1,2; &\quad & \widehat{u}_{3,j}=0,~~ j=1,2\\ & \widehat{u}_{i,3}=0,~~ i=0,1,2,3; &\quad & \widehat{u}_{i,0}= \sin(\pi\,i\,h_x)= \sin \frac{\pi\,i}{3}\,,~  i=0,1,2,3. \end{aligned}

С учетом краевого условия система может быть записана в следующей матричной форме


\begin{pmatrix}-4&1&1&0\\ 1&-4&0&1\\ 1&0&-4&1\\ 0&1&1&-4 \end{pmatrix}\!\cdot\! \begin{pmatrix}\widehat{u}_{1,1}\\ \widehat{u}_{1,2}\\ \widehat{u}_{2,1}\\ \widehat{u}_{2,2} \end{pmatrix}= \begin{pmatrix}-\sin \frac{\pi}{3}\\ 0\\-\sin \frac{2\pi}{3}\\0 \end{pmatrix}\!.

В результате имеем \widehat{u}_{1,1}= 0,\!325;~ \widehat{u}_{1,2}= 0,\!108;~ \widehat{u}_{2,1}= 0,\!325;~ \widehat{u}_{2,2}= 0,\!108.


Методика вычислений


1. Задать сетку D_h= \bigl\{(x_{i},y_{j})\bigr\}, где x_{i}= i\,h_x,~ y_{j}= j\,h_{y}, i=0,1,\ldots,I;~ j=0,1,\ldots,J; h_x,\,h_y — величины шагов по x и y;~ I\,h_x=L, J\,h_y= M,~ I и J — целые положительные числа; \varepsilon&gt;0 — желаемую точность. Вычислить


\begin{aligned}& \widehat{u}_{0,j}= \varphi_{0,j}, & & j=1,\ldots,J-1; &\quad & \widehat{u}_{I,j}= \varphi_{I,j}, & & j=1,\ldots,J-1;\\[2pt] & \widehat{u}_{i,J}= \varphi_{i,J}, & & i=0,1,\ldots,I; &\quad & \widehat{u}_{i,0}= \varphi_{i,0}, & & i=0,1,\ldots,I. \end{aligned}

Положить k=0; задать начальное приближение искомого решения во внутренних узлах: \widehat{u}_{i,j}^{(0)},~ i=1,\ldots,I-1; j=1,\ldots,J-1.


Разностная схема на сетке 3

2. Вычислить \widehat{u}_{i,j}^{(k+1)} одним из двух методов (рис. 8.13):


а) методом простых итераций: \widehat{u}_{i,j}^{(k+1)}= \frac{1}{4}\bigl(\widehat{u}_{i+1,j}^{(k)}+ \widehat{u}_{i-1,j}^{(k)}+ \widehat{u}_{i,j+1}^{(k)}+ \widehat{u}_{i,j-1}^{(k)}\bigr),\quad i=1,\ldots,I-1;\quad j=1,\ldots,J-1;

б) методом Зейделя:

Шаг 1. Положить j=1.

Шаг 2. Вычислить \widehat{u}_{i,j}^{(k+1)}= \frac{1}{4}\bigl(\widehat{u}_{i+1,j}^{(k)}+ \widehat{u}_{i-1,j}^{(k+1)}+ \widehat{u}_{i,j+1}^{(k)}+ \widehat{u}_{i,j-1}^{(k+1)}\bigr),\quad i=1,\ldots,I-1.

Шаг 3. Если j=J-1, процесс завершить. Иначе положить j=j+1 и перейти к шагу 2.


3. Если выполняется условие окончания \max_{i,j} \bigl|\widehat{u}_{i,j}^{(k+1)}-\widehat{u}_{i,j}^{(k)}\bigr| \leqslant \varepsilon, процесс завершить. Иначе положить k=k+1 и перейти к п.2.


Замечания


1. Метод Зейделя обеспечивает движение по узлам слева направо, начиная с первого слоя сетки по переменной x. Использование значений с (k+1)-й итерации, как правило, улучшает сходимость.


2. Вместо метода Зейделя можно использовать более эффективные релаксационные методы:


\widehat{u}_{i,j}^{(k+1)}= \omega\,\widehat{u}_{i,j}^{(k)}+ \frac{1-\omega}{4} \bigl(\widehat{u}_{i+1,j}^{(k)}+ \widehat{u}_{i-1,j}^{(k+1)}+ \widehat{u}_{i,j+1}^{(k)}+ \widehat{u}_{i,j-1}^{(k+1)}\bigr),

где \omega — параметр релаксации. Член \omega\, \widehat{u}_{i,j}^{(k)} характеризует "память" алгоритма о значении решения на предыдущей итерации. При \omega&gt;1 итерационный метод называется методом верхней релаксации, при 0&lt;\omega&lt;1 — методом нижней релаксации, при \omega= 0 совпадает с методом Зейделя. Существует оптимальное значение \omega^{\ast} параметра релаксации, зависящее от шагов сетки и способа получения \widehat{u}_{i,j}^{(k+1)}, при котором достигается максимальная скорость сходимости метода.


3. При решении задачи Неймана или третьей краевой задачи производные, входящие в краевые условия, следует также заменить их конечно-разностными аппроксимациями.


4. Если область \Omega, в которой решается задача, имеет неправильную форму, можно покрыть ее сеткой (рис. 8.14) и найти решение в ближайших к границе точках, а затем рассчитать значения искомой функции на границе, например, с помощью линейной интерполяции.


Область, покрытая сеткой

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

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

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


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

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