← Назад Далее →

3.2. Генерация разностной схемы

Для построения разностной схемы, соответствующей какому-либо шаблону, будем использовать интерполяционно-характеристический метод, описанный в п. 2.4. Каждому из узлов пространственно-временной сетки сопоставим значение сеточной функции ${\varphi }_i^n$. Набор узлов в шаблоне определяет, какие значения сеточной функции будут присутствовать в разностной схеме.

Рассмотрим сначала какую-либо явную схему. Например, схему «Крест», шаблон которой приведен на рисунке 30. На горизонтальной прямой $t=t_n$ в пространстве $(x,t)$ определим параметризацию $\xi =\frac{x-x_i}{x_{i+1}-x_i}$ и квадратичную функцию $f\left(\xi \right)=a_0+a_1\xi +a_2{\xi }^2$, коэффициенты которой пока не определены. Поскольку мы строим разностную схему для уравнения

$$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0,c=const,(3.1)$$

то решение этого уравнения будет постоянным вдоль характеристик - прямых $x+c\cdot t=const$.При выбранным шаблоне в разностной схеме участвуют четыре узла: $(x_i,t_{n-1})$, $(x_{i-1},t_n)$, $(x_{i+1},t_n)$ и $(x_i,t_{n+1})$, которым соответствуют значения сеточной функции $u_i^{n-1}$, $u_{j-1}^n$, $u_{i+1}^n$ и $u_i^{n+1}$. По значениям сеточной функции в трех точках можно определить коэффициенты функции $f(\xi )$. Два узла $(x_{i-1},t_n)$ и $(x_{i+1},t_n)$ уже лежат на прямой $t=t_n$ и определяют значения функции $f(\xi )$ в точках $\xi =-1$ и $\xi =1$. Еще один узел $(x_i,t_{n-1})$ лежит ниже прямой $t=t_n$. Выпустим из этого узла характеристику, которая пересечет прямую $t=t_n$ в точке $\xi =r$, где $r=\frac{c\cdot \tau }{h}$. Так как значение сеточной функции постоянно вдоль характеристики, то, тем самым, определено еще одно, третье, значение на прямой $t=t_n$(см. рис. 30).

Это приводит нас к системе из трех уравнений:

$$\{\begin{array}{c}f(-1)=a_0+a_1\cdot (-1)+a_2\cdot (-1{)}^2={\varphi }_{i-1}^n\\ f\left(1\right)=a_0+a_1\cdot (1)+a_2\cdot (1{)}^2={\varphi }_{i+1}^n\\ f\left(r\right)=a_0+a_1\cdot (r)+a_2\cdot (r{)}^2={\varphi }_i^{n-1}\end{array},(3.2)$$

разрешая которую, получаем коэффициенты $a_0$, $a_1<$ и $a_2$:

$$\{\begin{array}{l}{a}_0=\frac{{\varphi }_{i-1}^{n}r(r-1)+{\varphi }_{i+1}^{n}r(r+1)-2{\varphi }_i^{n-1}}{2(r^2-1)}\\ a_0=\frac{{\varphi }_{i+1}^n-{\varphi }_{i-1}^n}{2}\\ a_0=\frac{{\varphi }_{i-1}^n(r-1)+{\varphi }_{i+1}^n(r+1)-2{\varphi }_i^{n-1}}{2(r^2-1)}\end{array}.(3.3)$$

Рис. 30. Схема «Крест»

Таким образом, на прямой $t=t_n$ однозначно задана квадратичная функция $f(\xi )$. Учитывая, что характеристика, проходящая через четвертый узел шаблона $(x_i,t_{n+1})$ пересекает прямую $t=t_n$ в точке $\xi =-r$, разностная схема для уравнения (3.1) может быть записана в виде

$$f(-r)={\varphi }_i^{n+1}.(3.4)$$

Подставляя значения коэффициентов (3.3) в соотношение (3.4), находим:

$${\varphi }_i^{n-1}+r\cdot ({\varphi }_{i-1}^n-{\varphi }_{i+1}^n)={\varphi }_i^{n+1}(3.5)$$

или, в виде однородного уравнения:

$${\varphi }_i^{n-1}-{\varphi }_i^{n+1}+r\cdot ({\varphi }_{i-1}^n-{\varphi }_{i+1}^n)=0(3.6)$$

Соотношение (3.5) или эквивалентное ему (3.6) уже может рассматриваться как разностная схема, позволяющая получать значения сеточной функции на новом слое по времени $t_{n+1}$ по известным значениям сеточной функции на предыдущих слоях по времени $t_n$ и $t_{n-1}$.Такое представление схемы будем называть интерполяционным. Однако форма записи (3.6) не является окончательной и ниже будет показано, как эту форму записи можно улучшить.

Из сказанного выше следует, что разностную схему (3.6) определяют четыре уравнения (3.2) и (3.4) с тремя неизвестными $a_0$, $a_1$ и $a_2$, позволяющие, при исключении неизвестных, связать вместе одним соотношением значения сеточной функции в четырех узлах сетки. Это утверждение справедливо для любого шаблона, как для явных, так и для неявных схем. Тем самым, определен следующий алгоритм построения разностной схемы по заданному шаблону.

Для каждого узла $p_j,j=1...4$ ищется точка ${\xi }_j$ пересечения характеристики, проходящей через этот узел, с прямой $t=t_n$.

Решается система

$$\{\begin{array}{c}f({\xi }_1)=a_0+a_1{\xi }_1+a_2{\xi }_1^2={\varphi }_1\\ f({\xi }_2)=a_0+a_1{\xi }_2+a_2{\xi }_2^2={\varphi }_2\\ f({\xi }_3)=a_0+a_1{\xi }_3+a_2{\xi }_3^2={\varphi }_3\\ f({\xi }_4)=a_0+a_1{\xi }_4+a_2{\xi }_4^2={\varphi }_4\end{array}.$$

С помощью первых трех уравнений исключаются неизвестные $a_0$, $a_1$ и $a_2$, четвертое уравнение рассматривается как искомая разностная схема.

Фрагмент кода:

F[x_] := a0 + a1*x + a2*x^2; BuildScheme[m_] := Module[{eqs = {}, ic, jc, mi, x, ui, sch, cf},   For[ic = 3, ic >= 1, ic--,     For[jc = 1, jc <= 3, jc++,
mi = m[ [(ic - 1)*3 + jc] ];	(* узел сетки *)
If[mi != 0,	(* нашли точку нашего шаблона *)
x = (jc - 2) + (ic - 2) * r;	(* пересечение характеристики с опорной линией *)
ui = u[i + jc - 2, n - ic + 2];	(* значение сеточной функции в узле *)
eqs = Join[eqs, {ui == F[x]}],	(* добавляем к списку уравнение ui=F[x] *)
]	(* If[ip!=0 *)
]	(* For[jc=1 *)
];	(* For[ic=1 *)
sch = ui - Simplify[ui /. First @ Solve[conds, {a0, a1, a2, ui}]];   Return[sch]; ]; m = {0, 1, 0, 1, 0, 1, 0, 1, 0}; sch = BuildScheme[m]; Print[sch // TraditionalForm];

 

Out: -r u(i-1,n)+r u(i+1,n)-u(i,n-1)+u(i,n+1)

Определим порядок аппроксимации уравнения переноса (3.1) разностной схемой (3.6). В качестве точки разложения выберем точку, совпадающую с узлом $(x_i,t_n)$. Разложим входящие в схему (3.6) значения сеточной функции:

$${\varphi }_i^{n-1}={\varphi }_0-{\frac{\partial \varphi }{\partial t}|}_0\tau +{\frac{1}{2}\frac{{\partial }^2\varphi }{\partial t^2}|}_0{\tau }^2+o({\tau }^3);(3.7)$$ $${\varphi }_{i-1}^n={\varphi }_0-{\frac{\partial \varphi }{\partial x}|}_{0}h+{\frac{1}{2}\frac{{\partial }^2\varphi }{\partial x^2}|}_0{h}^2+o(h^3);(3.8)$$ $${\varphi }_{i+1}^n={\varphi }_0+{\frac{\partial \varphi }{\partial x}|}_{0}h+{\frac{1}{2}\frac{{\partial }^2\varphi }{\partial x^2}|}_0{h}^2+o(h^3);(3.9)$$ $${\varphi }_i^{n+1}={\varphi }_0+{\frac{\partial \varphi }{\partial t}|}_0\tau +{\frac{1}{2}\frac{{\partial }^2\varphi }{\partial t^2}|}_0{\tau }^2+o({\tau }^3),(3.10)$$

где индекс 0 означает значение в выбранной точке разложения. Подставим в (3.6) разложения (3.7)-(3.10). Слагаемые с четными степенями при $h$ и $\tau$сократятся. Для оставшихся нечетных степеней получим:

$$-2{\frac{\partial \varphi }{\partial t}|}_0\tau +r[-{2\frac{\partial \varphi }{\partial x}|}_{0}h+o(h^3)]+o({\tau }^3)=0.(3.11)$$

После вычитания (3.1) из (3.11) получим остаточный член

$$\begin{array}{l}-2{\frac{\partial \varphi }{\partial t}|}_0\tau +r[-2{\frac{\partial \varphi }{\partial x}|}_{0}h]+o(h^3,{\tau }^3)-{(\frac{\partial \varphi }{\partial t}+c\frac{\partial \varphi }{\partial x})}_0=\\ ={\frac{\partial \varphi }{\partial t}|}_0(-2\tau -1)+{\frac{\partial \varphi }{\partial x}|}_0(-2rh-c)+o(h^3,{\tau }^3)\end{array}$$

Отсюда видно, что разностная схема (3.6) аппроксимирует дифференциальное уравнение переноса (3.1) с нулевым порядком. С другой стороны, схема «крест» обладает вторым порядком аппроксимации.

Для устранения возникшего парадокса надо обратить внимание на то, что формально построенное соотношение (3.6) является однородным и может быть умножено на любое выражение. Заметим, что после подстановки в него разложений (3.7)-(3.10) мы получили выражение (3.11), в котором при производной по времени стоит коэффициент $-2\tau$, а в уравнении переноса - единица. Умножим (3.6) на $-\frac{1}{2\tau }$, учтем, что $r=\frac{c\cdot \tau }{h}$и получим традиционный вид разностной схемы:

$$\frac{{\varphi }_i^{n+1}-{\varphi }_i^{n-1}}{2\tau }+c\cdot \frac{{\varphi }_{i+1}^n-{\varphi }_{i-1}^n}{2h}=0(3.12)$$

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

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

Approx[sh_, center_, pow_] := Module[{sh1, sh2, ic, nc, cf},   Clear[i, n, x, t, h, tau];
ic = center[ [1] ]; nc = center[ [2] ];	(* центр аппроксимации *)
(* подставляем разложения в ряд Тейлора вместо сеточной функции: *)   sh1 = sh /. {     u[i, n] -> Series[u[x + (0 + (ic - i))*h, t + (0 + (nc - n))* tau], {h, 0, pow}, {tau, 0, pow}],     u[i + di_, n] -> Series[u[x + (di + (ic - i))*h, t + (0 + (nc - n))* tau], {h, 0, pow}, { tau, 0, pow}],     u[i, n + dn_] -> Series[u[x + (0 + (ic - i))*h, t + (dn + (nc - n))*tau], {h, 0, pow}, {tau, 0, pow}],     u[i + di_, n + dn_] -> Series[u[x + (di + (ic - i))*h, t + (dn + (nc - n))*tau], {h, 0, pow}, {tau, 0, pow}]   };
(* для вторых производных учитываем решение: *)   sh2 = sh1 /. {     D[(D[u[x, t], t] + c*D[u[x, t], x]), x] -> 0,     D[(D[u[x, t], t] + c*D[u[x, t], x]), t] -> 0   };
cfDuDt =Coefficient[sh2, D[u[x, t], t]];	(* коэффициент при производной по времени *)
cfDuDx = Coefficient[sh2, D[u[x, t], x]];	(* коэффициент при производной по пространству *)
Return[{cfDuDt, cfDuDx}]; ];

Модуль построения схемы BuildScheme, приведенный выше, следует закончить так:

...   sch = ui - Simplify[ui /. First @ Solve[conds, {a0, a1, a2, ui}]];   cf = Approx[sch, {i, n}, 2];   sch = Simplify[sch / cf[ [1] ] ];   Return[sch]; ];

← Назад Далее →