1252

Вычитая (7.48) из (7.49), получаем аппроксимацию

	da/dt |л+е да(an+1— an)/Atn				(7.50)
с погрешностью Е вида
Е =	{(1 - 0 )2 - 0 2 } (д/в/2) d*a/dt2|л+0 + ( Д е д {(1— 0)3 d3a/dt3|л+0>+
			-f-03d3a/d/s |п+0'}.		(7.51)
Отметим, что Е имеет порядок 0(At%)			при	0 = 1 /2 и 0(Д/„) в
противном случае. Умножив		равенство	(7.48)	на 1— 0, а	(7*49)
на 0 и сложив результаты, получим
	ал+0да(1 — 0)ал+ 0ал+1,				(7.52)
где	погрешность имеет порядок О (Д/3).
	Подставив в уравнение (7.47) аппроксимации (7.50) и (7.52) и
применив для fn+e то же представление,			что и для ал+е в		(7.52),
находим
{ ^	+ еК } аП+1 + { - д | + (1-	е) к } а л =	( 1 - 0 ) Г л+ 0 ^ .		(7.53)
Это соотношение в точности совпадает с равенством (7.46),					полу­
ченным методом конечных элементов.

Таким образом, с помощью подходящего выбора точки коллокации 0 внутри каждого элемента соотношение (7.46) можно прев­ ратить в одну из следующих известных конечно-разностных схем для уравнения (7.32).

1. Схема Эйлера (схема с разностью вперед). Эта		схема по­
лучается подстановкой в равенство (7.46) 0 =	0, что дает
(С/Дtn) an+1+ (— С/Д*я + К) ал =	V\	(7.54)

Если матрица С диагональна и, следовательно, элементы обрат­ ной матрицы С "1 могут быть получены непосредственно, то эта схема называется явной, так как ал+1 выражается в явном виде

через ал с помощью соотношения

ал+1 = Д ^С "1 {(С/Д/„ — К) а” + f"}.

(7.55)

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

схемы этого		типа	обычно сопровождаются		тем	недостатком,	что
используемый		шаг	по времени	Дtn не должен		превосходить	не­
которой	заданной		величины.
Матрица С будет диагональна				при конечно-разностной дискре­
тизации	пространственной области. Однако				если дискретизация
по пространству осуществляется				методом конечных элементов, то,

как правило, матрица С будет отлична от диагональной (см. ра­

венства (7.12)) и свойство явности этой схемы утрачивается. При

практических вычислениях данный недостаток обычно устраня­ ется посредством диагонализации матрицы С. Чтобы продемонст­

рировать этот процесс, рассмотрим задачу с одной пространст­

венной переменной и предположим, что при дискретизации по

пространству используются		одномерные линейные		элементы
с двумя узлами. Тогда	при	применении метода Галеркина при­
веденная матрица элемента се для			типичного элемента	е с номе­
рами узлов t и / задается равенством
	а(Щ)*		а
с е	а	а	(Щ)2_ dx.-	(7.56)

Если входящие в эту матрицу величины находить численным

интегрированием,	используя формулу	типа (5.10а) с	узлами на
концах элемента,	то
aNWj dx = (h'l2) [(а№,Щ)ужл, +		(aWf^ ) y3M,] =	0,	(7.57)
u*
так что матрица сг, а следовательно, и глобальная матрица				С диа-

гонализуются. Описание иных процедур диагонализации и под­

робное их обсуждение можно найти в книге [2].
2. Схема		Кранка—Николсона (схема с центральной разностью).
При 0 = 1 / 2		равенство (7.46)			принимает вид
( k + 7 К ) а ” ‘ + (				'	, + Т К ) а " = Т Р '+ Т » * * ■ •						( 7 - 5 8 )
Говорят, что эта схема является неявной, так как									для	нахожде­
ния a'i+1 требуется			решить	систему		уравнений, матрица				которой
отлична	от диагональной.
3. Схема с разностью назад. Эта схема получается подста­
новкой	в равенство		(7.46)	0 = 1, что			приводит к		неявной		схеме
		(С/Д/я + К) а”+1— (С/Дtn) а” =						in+1.			(7.59)
Очевидной		альтернативой		поточечной			коллокации		является		при­
менение на каждом			элементе метода			типа Галеркина, т. е. исполь­
зование на элементе п равенства						W.n = N% или Wn= /V"+1. Этот
способ	введения		весовых		функций		дает	в точности		частный
случай	равенства		(7.46).	Читатель			может	показать,		что	при
=		получающаяся схема идентична схеме получаемой при­
менением поточечной коллокации при Т = 2/3								на каждом элементе.

7.5.3.Многослойные схемы для уравнений первого порядка

Если отказаться от ограничения (7.37) и рассматривать весо­

вые функции Wn, отличные от нуля более чем на одном элемен­ те, то можно получить схемы с любым числом слоев. Чтобы продемонстрировать этот процесс, вернемся к равенству (7.36) и

непосредственно воспользуемся стандартным подходом Галеркина,

т. е. положим
Так как		(7.60)
yv„ = 0 при	и	(7.61)
то уравнение метода взвешенных невязок принимает вид
in+1
J (Cda/ Я + К а — Uf))Nndt = 0 ,	n = 1 , 2 , 3 , . . .	(7.62)
/п-1

Применив используемый для а вид интерполяции также и для f

(см. равенство (7.43)),

получим

<и+.

{С (a"-1 dNn_Jdt +

a" dNjdt) +

dt +

j

К (

a

+

a"iVn)}

in- i

in+1

an+1 dNn+1ldt) +

a"+Wn+1)} Nn dt =

+

J

{C (a" diV„/df +

К (a"W„ +

in

=

j

a n- 1^ n- 1 + f ^ n) ^ nd / +

jf

(f»^„ +

f » + ^ B+l) ^ Bdf.

(7.63)

-1

in

Используя соотношения

(7.35)

и проводя

интегрирование,

полу­

чаем трехслойную схему

fп + 1

2Vl .

f”- 1

. . * . + ! K

. . +

( - 5

§ j + * ) . -

’ (7.64)

6

+

3

6

где

предполагается,

что

Atn_1 = htn = kt1 а С и

К — постоянные

матрицы. Вычисления организуются с помощью повторного использования такой схемы при условии, что известны два на­ чальных значения а0 и а1. Так как а0 задается с помощью на­

чальных данных, то требуется каким-либо образом вычислить а1.

Это можно сделать, например, использовав на исходном шаге одну из описанных выше двухслойных схем. Дальнейшие сведе­ ния о других схемах этого типа приводятся в [2].

Пример 7.4. Применим описанные выше схемы_для интегриро­

вания простейшего скалярного уравнения da/dt +					ка = 0 (к— поло­
жительная		постоянная) с		начальным условием	а = 1 в момент
времени	t =	0.
В этом случае общее уравнение (7.46) для двухслойной схемы
упрощается,		принимая	вид
( t k	+ e ) “ ” t , + ( ~		t k	+ ( 1 _ 9 ) ) “ ’ " 0 ’ п - ° ' ' • 2>

причем согласно начальному условию а(/ = 0) = а° = 1.

Рис. 7.5. Результат применения некоторых двухслойных схем к уравнению daldt-\- + Я а = 0 при ЯД/„=0.5.

Рис. 7.6. Результат применения некоторых двухслойных схем к уравнению da/dt+ -+-Яа=0 при Я,Д/„=2.2. Обозначения те же, что на рис 7.5.

На рис.

7.5 и

7.6

поведение точного

решения сравнивается

с поведением

решения,

полученного по этой

схеме

при 0 =

0 (раз­

ность вперед),

0 =

1/2

(схема

Кранка—Николсона), 0 = 2/3

(метод

Галеркина)

и 0 = 1

(разность

назад),

для двух различных

шагов

по времени. Согласно

рис. 7.5, результаты,

полученные

по всем

схемам,

дают

правильную

картину поведения решения, когда шаг

по времени

удовлетворяет

соотношению

^Д/п = 0.5,

причем наи­

более точные

результаты

получаются

при 0 = 1 /2 , т.

е. по схеме

Кранка—Николсона.

Этого

и следовало ожидать,

так

как

не­

трудно

проверить,

что

погрешность

аппроксимации da/dt |п+е в

этой схеме имеет порядок О(Д^), а

для других

значений 0

по­

рядок 0(Atn) (см. равенство (7.51)).

Согласно

рис.

7.6,

при возрастании

шага

по

времени

до

зна­

чения,

определяемого

соотношением

А,Д^ = 2.2,

поведение реше­

ния для некоторых

схем

заметно изменяется. Решение,

получае­

мое по схеме с разностью вперед, дает колебания с нарастающей

амплитудой. В этом случае говорят, что для данного значения

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

Николсона, хотя и остающаяся устойчивой, утрачивает свою вы­

сокую точность и обнаруживает склонность к колебаниям, отсут­ ствующим в точном решении. Метод Галеркина и схема с раз­

ностью назад продолжают давать профили правильной формы

с заметным снижением точности по сравнению с получаемой при

меньшем значении шага по времени, как и следовало ожидать согласно (7.51). Обсуждение свойств устойчивости решений схем

интегрирования по времени и возможности возникновения коле­ баний будет продолжено в следующем разделе.

Трехслойная схема (7.64) для этой задачи принимает вид

и опять согласно начальным данным а° = 1.

Чтобы вычислить значение а1 и тем самым получить возмож­

ность рекуррентно применять трехслойную схему, можно вос­

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

нением трехслойной		схемы, на исходном шаге которой исполь­
зуется двухслойная		схема типа разности вперед,	с шагом по
времени,	определяемым соотношением XДtn= 0.1,		сравниваются
с точным	решением. Удается получить хорошую		точность при­

ближения, хотя можно заметить, что найденные числовые значе­ ния колеблются около точного решения.

It	Т о ч н о е	Т р е х с л о й н а я
	р е ш е н и е	с х е м а
0	1.0	1.0
0.1	0.9048	0.9
0.2	0.8187	0.8194
0.3	0.7408	0.7362
0.4	0.6703	0.6715
0.5	0.6065	0.6021

7.5.4. Характеристики устойчивости двухслойных схем

Общее уравнение первого порядка (7.32) с нулевой правой частью f может быть решено численно с помощью двухслойной схемы (7.53), записываемой в виде

(1/Дг„ +		е с ^ К ) an+1 =				(1 — 0) С ^К ) а".					(7.65)
Если Х1( Х2,			— собственные значения,					а	ах,	а2,	ам—
соответствующие		собственные			векторы		задачи		на	собственные
значения (см.	разд.		7.4.2)	Кос	ХСсс =	0,					(7.66)
то в общем случае мы имеем
			м			м
	ат =		т2—1Ута а >		ап+1=	т2= I Ут+'ат‘,					(7.67)
подставив эти	выражения			в (7.65), получим
			Ут ~	\/Мп- М т			Ут-				(7.68)
Как было отмечено выше, в случае						положительно				определенных
матриц К и С все			собственные значения Xlf					h2f		Км задачи
(7.66) вещественны			и неотрицательны.
Следовательно, точное решение затухает со временем и, со­
гласно (7.67), то же самое поведение						будет характерно					для чис­
ленного решения,		если
т. е. если			I < I f/m I.		tn =l,	2,		М,			(7.69)
- 1	<	- - иА| п-0хе^			И< '■						(7.70)

Зто условие должно быть выполнено для устойчивости схем интегрирования по времени. Кроме того, схема будет свободна от

колебаний, если каждая входящая в виде множителя в у\1п мода

имеет один и тот же знак на всех временных слоях /г. В силу

равенства (7.68) это будет выполняться, если

Ут	1/А<в + ( 1 —	8) Хт	„	т = 1,2,	М.	(7.71)
Упт	l / A t n - Q	\ m	> и ’

Следовательно, двухслойная по времени схема (7.65) будет устой­

чива и свободна от колебаний, если

	0 ^	— 0) Хт^ 1	т— 1, 2,	М.	(7.72)
	^	1 / д < „ - е \ и ^

Условие (7.70) или (7.71) при заданном значении 0 обычно будет вводиться как ограничение на максимальный размер шага по вре­

мени, который может быть использован в той или иной конк­

ретной задаче.

Пример 7.5. Рассмотрим свойства устойчивости двухслойных

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

рено решение одного уравнения da/dt-\-Xa = 0 с использованием схемы вида

(1/Д/„ +		Х0) ап+1 +	(-1/Д/„ + X (1-0)) а"=0.
Для данного	случая одного		уравнения	решение	системы (7.66^
тривиально и	мы	имеем Х =	— X. Тогда	условие	устойчивости
(7.70) сводится к двойному неравенству
		1/Д^—М Ь-0)		j
		^	1/д * „ + ).е

Второе из этих неравенств автоматически выполняется, а для выполнения первого требуется, чтобы

(1 — 2 0 ) Х д / п < 2 .

Это условие всегда выполняется для 0 ^ 1/2 каково бы ни было

значение Д О таких схемах говорят, что они являются безус­ ловно устойчивыми. При О < 0 < 1 / 2 устойчивость является услов­

ной, и требуется, чтобы шаг по времени Д и сп о л ь з у е м ы й на

каждом этапе процесса, удовлетворял ограничению

ЯД / „ < 2Д 1-20).

Вданном случае схема, свободная от колебаний, может быть

получена при условии, что выполняется соотношение (7.72), т. е.

если

1/А/п-Х(1_-0) р 1/Д^л + ^б

а это эквивалентно требованию

( 1 - в ) К М п < 1.

Данное условие всегда выполняется			для схемы с разностью назад
0 = 1	каково бы ни было значение шага по времени Дtn. Однако
для	всех других	значений 0 при достаточно больших значениях
шага	по времени	будут возникать	колебания. Из приведенного

неравенства следует, что использованный при анализе примера 7.4

шаг по времени, определяемый соотношением ЯД/П= 2.2, таков,

что обеспечено возникновение колебаний для схемы Кранка —

Николсона (0=1/2), но при этом метод Галеркина (0 = 2/3) будет

свободен от них (см. рис. 7.6). Таким образом, при использовании

достаточно больших значений шага по времени математически

более точная схема Кранка — Николсона на практике может давать

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

колебаний в плоскости (0, ЯД/„) для	рассматриваемых			в данном
случае двухслойных схем.
В разд. 7.4.3 было показано, как	с	помощью разложения по
модам перейти от векторного уравнения
С da/d/+ Ка =	О			(7.73)
к системе независимых уравнений вида (7.30):
dylldt— 'klyl^ 0, k=\ ,	2,	3,	М.	(7.74)

Следовательно, можно получить информацию об ожидаемом пове­ дении решения, получаемого применением к уравнению (7.73) двухслойной схемы, рассматривая применение этой схемы к каж­ дому уравнению вида (7.74). Это в свою очередь можно сделать, повторив проведенные в примере 7.5 рассуждения для каждого из этих независимых уравнений. Можно ожидать, что решение

будет свободно от колебаний, если

Д * „ < 1 / [ ( 1 - 0 ) 1 М ] ,

и будет устойчивым, если

(1 — 20) | Xt | Д/„ < 2

(7.76)

для любого значения К,. Если используется безусловно устойчи­ вая схема (т. е. выбирается такое 0, что неравенство (7.76) выпол­ няется автоматически) и некоторый шаг по времени АГп, только частично удовлетворяющий неравенству (7.75) {например, шаг по времени Д/л* может быть таким, что

Д / ; < 1 /[(1 - 0 )|^ П , 1=1, 2 , . . . , р ,

(7.77)

НО

Д / ; > 1 /[(1 -0 )|Х ,|], / = р + 1 , р + 2.........

М),

(7.78)

Рис. 7.7. Области неустойчивости и колебаний для двухслойных схем, применяе­

мых	к уравнению		da/dt-\-ka= 0.
то	моды для	малых значений \Kt \ будут свободны от колебаний,
но	для мод	с	большими значениями |XZ| колебания будут воз­

никать. Очевидным выходом из положения было бы уменьшение

шага по времени	до такого значения,		чтобы неравенство (7.75)
выполнялось для	всех	однако для	этого может потребоваться

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

мод. Обычно такие процедуры состоят в каком-либо осреднении

вычисленных значений последовательно на каждом шаге по вре­

мени.

Условно устойчивые схемы всегда будут приводить к неустой­

чивости, если используемый шаг по времени Д/„ не удовлетворяет неравенству (7.76), в котором за 'к1 взято максимальное по абсо­

лютной величине собственное значение задачи (7.66). Это сущест­ венно ограничивает использование явного метода (7.55). Однако

необходимость использования при применении этого метода доста­

точно большого числа малых шагов по времени часто компенси­

руется ненадобностью обращения матрицы, что существенно

экономит вычислительные усилия.

Согласно проведенному в этом разделе анализу устойчивости

процесс получения оценки для максимального шага по времени,

который может быть использован в условно устойчивой схеме,

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

значения (7.66):

Ксс+^Са = 0.

(7.79)

Однако можно показать [3], что если дискретизация по прост­

ранству осуществляется методом конечных элементов, то наиболь­

шее глобальное собственное значение всегда должно быть меньше

наибольшего локального собственного значения элемента. Эти

собственные значения элемента нетрудно получить из уравнения (7.79), заменив в нем глобальные матрицы приведенными матри­

цами элемента. Тогда можно провести оценку максимального

допустимого шага по времени. Следует заметить, однако, что

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

Пример 7.6. Используем двухслойную схему типа разности

вперед (т. е. с 0 = 0) для решения уравнения Cda/d/ + Ka = f, получающегося при применении метода конечных элементов по

пространственной переменной к задаче из примера 7.2. Тогда ненулевые компоненты типичных приведенных матриц элементов се и ке имеют вид

1/3	1/6“	1	— Г
сe = he 1/6	1/3J ’ ке	he — 1	1 •

В данном случае задача на собственные значения для элемента записывается как

det (k* + Xc*) = 0.

Подставляя сюда приведенные выше выражения для матриц эле­

ментов, получаем уравнение

Г l + M/i*)2/3 — 1 + Х(А*)в/6]

[ — 1 + к ( И 2/6 1 + х ( И 2/з] ^ ° ’

решение которого дает

*1 = 0, Х2 = -12/(А ‘)».

Таким образом, условие устойчивости (7.76) выполняется, если

2) В случае равномерной сетки эта оценка иногда не допускает сущест­ венного улучшения.— Прим, ред.