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

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

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

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

$$\frac{\partial u}{\partial t} +c\frac{\partial u}{\partial x} =0,\quad 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({\rm\xi })\). Два узла \((x_{i-1} ,t_{n} )\) и \((x_{i+1} ,t_{n} )\) уже лежат на прямой \(t = t_n\) и определяют значения функции \(f({\rm\xi })\) в точках \({\rm\xi }=-1\) и \({\rm\xi }=1\). Еще один узел \((x_{i} ,t_{n-1} )\) лежит ниже прямой \(t = t_n\). Выпустим из этого узла характеристику, которая пересечет прямую \(t = t_n\) в точке \({\rm\xi }=r\), где \(r=\frac{c\cdot {\rm\tau }}{h}\). Так как значение сеточной функции постоянно вдоль характеристики, то, тем самым, определено еще одно, третье, значение на прямой \(t = t_n\)(см. рис. 30).

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

$$\left\{\begin{array}{c} {f\left(-1\right)=a_{0} +a_{1} \cdot (-1)+a_{2} \cdot (-1)^{2} ={\rm\phi }_{i-1}^{n} } \\{f\left(1\right)=a_{0} +a_{1} \cdot (1)+a_{2} \cdot (1)^{2} ={\rm\phi }_{i+1}^{n} } \\{f\left(r\right)=a_{0} +a_{1} \cdot (r)+a_{2} \cdot (r)^{2} ={\rm\phi }_{i}^{n-1} } \end{array}\right. , (3.2)$$

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

$$\left\{\begin{array}{l} {a_{0} =\frac{{\rm\phi }_{i-1}^{n} r\left(r-1\right)+{\rm\phi }_{i+1}^{n} r\left(r+1\right)-2{\rm\phi }_{i}^{n-1} }{2(r^{2} -1)} } \\{a_{0} =\frac{{\rm\phi }_{i+1}^{n} -{\rm\phi }_{i-1}^{n} }{2} } \\{a_{0} =\frac{{\rm\phi }_{i-1}^{n} \left(r-1\right)+{\rm\phi }_{i+1}^{n} \left(r+1\right)-2{\rm\phi }_{i}^{n-1} }{2(r^{2} -1)} } \end{array}\right. . (3.3)$$

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

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

$$f(-r)={\rm\phi }_{i}^{n+1} . (3.4)$$

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

$${\rm\phi }_{i}^{n-1} +r\cdot \left({\rm\phi }_{i-1}^{n} -{\rm\phi }_{i+1}^{n} \right)={\rm\phi }_{i}^{n+1} (3.5)$$

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

$${\rm\phi }_{i}^{n-1} -{\rm\phi }_{i}^{n+1} +r\cdot \left({\rm\phi }_{i-1}^{n} -{\rm\phi }_{i+1}^{n} \right)=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\) ищется точка \({\rm\xi }_{j}\) пересечения характеристики, проходящей через этот узел, с прямой \(t = t_n\).

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

$$\left\{\begin{array}{c} {f({\rm\xi }_{1} )=a_{0} +a_{1} {\rm\xi }_{1} +a_{2} {\rm\xi }_{1} ^{2} ={\rm\phi }_{1} } \\{f({\rm\xi }_{2} )=a_{0} +a_{1} {\rm\xi }_{2} +a_{2} {\rm\xi }_{2} ^{2} ={\rm\phi }_{2} } \\{f({\rm\xi }_{3} )=a_{0} +a_{1} {\rm\xi }_{3} +a_{2} {\rm\xi }_{3} ^{2} ={\rm\phi }_{3} } \\{f({\rm\xi }_{4} )=a_{0} +a_{1} {\rm\xi }_{4} +a_{2} {\rm\xi }_{4} ^{2} ={\rm\phi }_{4} } \end{array}\right. .$$

С помощью первых трех уравнений исключаются неизвестные \(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) значения сеточной функции:

$${\rm\phi }_{i}^{n-1} ={\rm\phi }_{0} -\left. \frac{\partial {\rm\phi }}{\partial t} \right|_{0} \tau +\left. \frac{1}{2} \frac{\partial ^{2} {\rm\phi }}{\partial t^{2} } \right|_{0} {\rm\tau }^{2} +o({\rm\tau }^{3} ); (3.7)$$ $${\rm\phi }_{i-1}^{n} ={\rm\phi }_{0} -\left. \frac{\partial {\rm\phi }}{\partial x} \right|_{0} h+\left. \frac{1}{2} \frac{\partial ^{2} {\rm\phi }}{\partial x^{2} } \right|_{0} h^{2} +o(h^{3} ); (3.8)$$ $${\rm\phi }_{i+1}^{n} ={\rm\phi }_{0} +\left. \frac{\partial {\rm\phi }}{\partial x} \right|_{0} h+\left. \frac{1}{2} \frac{\partial ^{2} {\rm\phi }}{\partial x^{2} } \right|_{0} h^{2} +o(h^{3} ); (3.9)$$ $${\rm\phi }_{i}^{n+1} ={\rm\phi }_{0} +\left. \frac{\partial {\rm\phi }}{\partial t} \right|_{0} {\rm\tau }+\left. \frac{1}{2} \frac{\partial ^{2} {\rm\phi }}{\partial t^{2} } \right|_{0} {\rm\tau }^{2} +o({\rm\tau }^{3} ), (3.10)$$

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

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

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

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

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

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

$$\frac{{\rm\phi }_{i}^{n+1} -{\rm\phi }_{i}^{n-1} }{2{\rm\tau }} +c\cdot \frac{{\rm\phi }_{i+1}^{n} -{\rm\phi }_{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]; ];

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