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

− f (ti , x j )− ft′(ti , x j )τ2 + fx′(ti , x j )h2 + O(τ2 ,h2 )=

+ [fx′(ti , x j )− vu′xx′ (ti , x j )]h2 + O(τ2 ,h2 ).

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

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

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

В силу того, что αmax =1τ, для устойчивости по начальным данным требуется выполнение условия

1τ ≥ 1 τ − v h + vh ,

1τ − v h ≤1 τ − v h ,

−1τ + vh ≤1τ − vh ≤1τ − vh ,

2vh ≤ 2τ.

Разностная схема условно устойчива при выполнении условия Куранта1, ограничивающего шаги интегрирования:

Устойчивость по правой части имеет место при выполнении условия

1τ ≥ ωτ,

что справедливо, например, при ω ≤1. Из выражения (11.4) следует разрешимость и единственность решения разностного уравнения (11.3),

1 Курант Рихард [8.1.1888 – 27.1.1972] – математик, учился в университетах Бреслау и Цюриха. Профессор Геттингенского университета с 1920 по 1933 годы, профессор НьюЙоркского университета с 1934 года. Основные научные труды выполнены в теории конформных отображений и для краевых задач математической физики. В 1966 году был избран иностранным членом Академии наук СССР.

251

u j = (1− ντh)u j + ντhu j−1 + f τ.

Согласно теореме 10.4 численное решение задачи (11.1)–(11.2) на разностной схеме (11.3) сходится к точному решению при выполнении соотношения (11.5) и уменьшении шагов по времени и координате, то есть при h, τ → 0.

Вторая разностная схема, соответствующая шаблону, показанному на

и также имеет погрешность аппроксимации O(h, τ). Условие устойчивости по начальным данным в этом случае записывается иначе:

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

Разностная схема для шаблона, изображенного на рис. 11.1, в,

и имеет также первый порядок аппроксимации.

Шаблон, показанный на рис. 11.1, д, приводит к разностному аналогу дифференциального уравнения (11.1)

имеющему второй порядок погрешности аппроксимации.

Формально разностные схемы, получаемые для шаблонов, изображенных на рис. 11.1, а, и 11.1, г, являются явными, а на рис. 11.1, б, 11.1, в и 11.1, д – неявными. В то же время в трех последних случаях можно явно вычислить зна-

чение u j . Действительно, крайнее левое значение u0 =U 0 (t ) известно из граничного условия (11.2). Следовательно, можно определить u1, рассматривая

разностную схему как уравнение с одним неизвестным. Затем можно найти u2 , пользуясь разностным уравнением, записанным для следующего узла, и так далее для всех остальных узлов сетки.

Схема (11.9) имеет второй порядок аппроксимации и при не слишком больших шагах h и τ дает более точное решение. Однако на быстропеременных решениях лучшиерезультатыпоточностиполучаютсянасхемах(11.3), (11.6) и(11.8).

252

Разностная схема (11.8) безусловно устойчива, что имеет преимущество при проведении вычислительного эксперимента. В то же время схемы (11.3) и (11.6) являются более точными.

11.1.2. Явно-неявная схема

Пусть шаги интегрирования h и τ являются переменными, а скорость движения частицы v(t, x) зависит от времени t и координаты x. В этом случае при использовании разностных схем (11.3) и (11.6) для определения решения на очередном шаге интегрирования необходимо выполнение условия Куранта для устойчивости процесса вычислений. Целесообразно строить вычислительный процесс следующим образом: если выполняется условие (11.5), то очередное значение uˆ j определять из разностной схемы (11.3); при справедливости усло-

вия (11.7) вычисления выполнять по формуле (11.6). Это гарантирует устойчивость явно-неявной схемы; следует отметить ее более высокую точность по сравнению со схемой (11.8).

11.1.3. Схема для двумерного уравнения переноса

Рассматривается задача о переносе частицы потоком с двумя пространственными переменными, x и y,

253

Для оценки порядка погрешности аппроксимации дифференциального уравнения (11.10) разностной схемой (11.12) используются разложения решения и правой части уравнения в ряды Тейлора:

u(ti−1, x j , yk )= u(ti , x j , yk )− ut′(ti , x j , yk )τ + utt′′(ti , x j , yk )τ2 2 + O(τ3 ),

u(t , x − , y )= u(t , x , y )− u′ (t , x , y )h + u′′ (t , x , y )h2 2 + O(h3 ),

i j 1 k i j k x i j k x xx i j k x

u(t , x , y − )= u(t , x , y )− u′ (t , x , y )h + u′′ (t , x , y )h2 2 + O(h3 );

i j k 1 i j k y i j k y yy i j k y

f(ti−1 2 , x j−1 2 , yk−1 2 )=

=f (ti , x j , yk )− ft′(ti , x j , yk )τ2 − fx′(ti , x j , yk )hx 2 − f y′(ti , x j , yk )hy 2 + O(τ2 ,hx2 ,hy2 ).

Эти разложения подставляются в разностный аналог (11.12):

ψ jk = ut′(ti , x j , yk )−utt′′(ti , x j , yk )τ2 +u′x (ti , x j , yk )vx −u′xx′ (ti , x j , yk )hx 2 +

+ fx′(ti , x j , yk )hx 2 + f y′(ti , x j , yk )hy 2 + O(τ2 ,hx2 ,hy2 )=

= [ut′(ti , x j , yk )+ u′x (ti , x j , yk )vx + u′y (ti , x j , yk )vy − f (ti , x j , yk )]+

+ [ft′(ti , x j , yk )− utt′′(ti , x j , yk )]τ2 + [fx′(ti , x j , yk )− u′xx′ (ti , x j , yk )]hx 2 +

+ [f y′(ti , x j , yk )− u′yy′ (ti , x j , yk )]hy 2 = O(τ,hx ,hy ).

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

u jk (1 τ + νx hx + νy hy ) − u j−1kνx hx − u jk−1νy hy = u jk τ + f .

В этом выражении

αmax =1τ + vx hx + vy hy.

Всоответствии с принципом максимума условие устойчивости по начальным данным требует выполнения неравенства

1τ + vx hx + vy hy ≥ vx hx + vy hy +1τ,

254

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

(1τ + vx hx + vy hy )− (vx hx + vy hy )=1τ ≥ ωτ при ω ≤1.

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

11.2. Уравнения параболического типа

Линейное дифференциальное уравнение

описывает процесс нестационарной теплопроводности в тонком однородном стержне либо диффузии. В первом случае коэффициент η играет роль коэффициента температуропроводности, во втором случае – коэффициента диффузии.

Для оценки порядка погрешности аппроксимации дифференциального уравнения (11.13) разностным аналогом (11.15) точное решение разлагается в ряды Тейлора1 возле точки (ti , x j ),

1 Ряд Тейлора для функций нескольких переменных (кратный ряд):

255

f (ti+12 , x j )= f (ti , x j )+ ft′(ti , x j )τ2 + ftt′′(ti , x j )τ2 8 + O(τ3 ), u(ti+1 , x j+1 )= u(ti , x j )+ut′(ti , x j )τ+u′x (ti , x j )h +

+[utt′′(ti , x j )τ2 + 2utx′′ (ti , x j )τh +u′xx′ (ti , x j )h2 ]2 +

+[uttttiv (ti , x j )τ4 + 4utttxiv (ti , x j )τ3h + 6uttxxiv (ti , x j )τ2 h2 + + 4utxxxiv (ti , x j )τh3 +uxxxxiv (ti , x j )h4 ]24 +O(τ5 , h5 ),

u(ti+1 , x j−1 )= u(ti , x j )+ut′(ti , x j )τ−u′x (ti , x j )h + + [utt′′(ti , x j )τ2 − 2utx′′ (ti , x j )τh + u′xx′ (ti , x j )h2 ]2 +

+ uttttiv (ti , xj )τ4 − 4utttxiv (ti , xj )τ3h + 6uttxxiv (ti , xj )τ2h2 −

−4utxxxiv (ti , xj )τh3 + uxxxxiv (ti , xj )h4 − / 24 + O(τ5 ,h5 ).

В результате подстановки этих разложений в формулу (11.15) получается оценка погрешности

Сучетом уравнения (11.13) это выражение дает оценку

ψij = [utt′′(ti , x j )− 2ησutxx′′′ (ti , x j )− ft′(ti , x j )]τ2 +

256

+[4uttt′′′(ti , x j )−12ησuttxxiv (ti , x j )−3 ftt′′(ti , x j )]τ2 24 −

− uxxxxiv (ti , x j )ηh2 12 + O(τ3 ,h3 )= O(τ,h2 ),

то есть имеет место погрешность аппроксимации первого порядка по шагу τ и второго порядка по шагу h. Если продифференцировать по переменной t вы-

ражение (11.13),

utt′′ = ηutxx′′′ + ft′(t, x),

и положить весовой коэффициент σ =12 , то погрешность аппроксимации можно понизить до второго порядка по обоим шагам интегрирования:

свесовым коэффициентом σ =12 носит название схемы Крэнка–Николсона. Для оценкиустойчивостиразностнойсхемывыражение(11.15) преобразуетсяквиду

u j−1 (−ησh2 ) + u j (1 τ + 2ησ h2 ) + u j+1 (−ησh2 ) =

Выбирается коэффициент

αmax =1τ + 2ησh2 .

Согласно принципу максимума для устойчивости решения по начальным данным должно выполняться условие

1 τ+ 2ησ h2 ≥ 2ησ h2 + 2η(1−σ) h2 + 1 τ−2η(1−σ)h2 ,

1τ − 2η(1 − σ) h2 ≥ 1 τ − 2η(1 − σ)h2 ,

2η(1 − σ)h2 −1 τ ≤1 τ − 2η(1 − σ)h2 ≤1τ − 2η(1 − σ)h2 ,

2η(1 − σ)h2 ≤1τ,

τ ≤ h2 2η(1 − σ).

Устойчивость по правой части имеет место в случае неравенства

257

1τ + 2ησh2 − 2ησh2 ≥ ωτ, ω> 0 ,

которое выполняется, например, при ω<1. Рассмотренная схема имеет первый (или второй при σ = 1/2) порядок аппроксимации и является условно устойчивой по начальным данным и по правой части. Это означает, что семейство рассмотренных схем с весами дает последовательности численных решений, сходящихся к точному решению при уменьшении шагов τ и h и выполнении указанного соотношения между шагами τ и h. Для схемы Крэнка–Николсона условие устойчивости выполняется, если

τ≤ h2 η.

11.2.2.Трехслойная схема Ричардсона

Для уравнения (11.13) разностный аналог с использованием пятиточечного шаблона (рис. 11.4) записывается в форме

[(ρk )m+1 − (ρk )m−1 ]eikx j 2τ = (ρk )m η[eik (x j −h) − 2eikx j + eik (x j +h)]h2,

ρk2 −ρk 2τη(e−ikh − 2 + eikh )h2 −1= 0.

Использование формулы Эйлера приводит к выражению

ρ2k + ρk 8τηsin2 (kh2)h2 −1 = 0.

Корни этого уравнения вещественные, различные,

ρk = 4τηsin2 (kh2)h2 ± 16τ2η2 sin4 (kh2) h4 +1 .

258

Поскольку произведение корней

(ρk )1 (ρk )2 =1,

можно сделать вывод о том, что один из них превышает единицу. Но согласно теореме 10.5 это означает неустойчивость разностной схемы (11.16).

11.2.3. Схема Дюфорта и Франкела

На том же шаблоне (см. рис. 11.4) можно построить другой разностный аналог уравнения (11.13):

зависит от соотношения шагов интегрирования τh и мала в случае

lim τ h = 0.

τ,h→0

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

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

259