Здесь х{+ ~ йух|+ -f (1 - w)xi (аналогично определяется г/+°).

Если аналитическое дифференцирование fag оказывается сложным, используются процедуры экстраполяции, при которых значение х'. * для вычисления Ц^ получают экстраполяцией, исходя из значений на нескольких предыдущих слоях (с верхними индексами /, /- 1, ...). Например, полагают

х{+ {1+w)xl - wxb\ (6.4.11)

jc+ ~ [1 -f 1,5йу -f 0,5йу2] xi - (йу2 -f 2w) xr + 0,5(w- + w) x{-.

(6.4.12)

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

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

> Здесь имеется в виду, что f(x, у) (или g(x, у)) содержит слагаемое х^, и обсуждается, какое слагаемое можно ввести в определение величины /+ . - Прим. ред.

2) Разумно сначала определить по (6.4.7) и потом оставить тольКс* первые члены тейлоровского разложения. По существу об этом уже шли речь (см. абзац после формулы (6.4.7)). - Яриж ред.

) Здесь имеется в виду, что (по определению) = f (jcJ-+° , yi ) Лрим. ред.

(л;+)2 МОЖНО заменить) произведением (xi+i xi); точно так же {x{+f{x{f х1.+К Приау=1/2 член {xi+f можно заменить

выражением [{х{ + х{+)/2у ~ [{xl.f -f 2л;,Ц+ + xixi+]. Довольно часто используется еще более простая замена

fT-f[. (6.4.9)

Если функции / и легко продифференцировать аналитически, то можно воспользоваться формулой Тейлора, т. е. (например, для fi+ ) положить) (ср. с (6.4.7))

6.4.2. Автоматическое изменение шага по времени

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

a) Сравнение значений решения на новом слое, полученных с шагом т (х'+с^), со значениями, полученными с помощью-двух шагов длиной т/2 {х'+-Ч^>). В этом случае оценкой погрешности с целью регулирования шага может служить величина

- II л;/+ W - л;/+ II (6.4.13)

(аналогичная формула записывается для Еу). Для более тонкой оценки Е можно также использовать экстраполяцию по Ричардсону [6.10]. Если действовать так, то для одного удачного и контролируемого шага т нам потребуется сделать три пробных шага. По этой причине указанный подход используется сравнительно редко и в большинстве случаев после нескольких постоянных (т. е. не регулируемых) шагов по времени.

b) Оценка различия между экстраполированным значением найденным с помощью формул (6.4.11) или (6.4.12), и

вычисленным значением

Ех = тах xi+ - wx{+ - (\ - W) х1\. (6.4.14)

Аналогично определяется Еу. Такую же оценку можно получить и для метода предиктор-корректор .

c) В случае использования какой-то схемы линеаризации (6.4.8) мы можем положить

= max I / {wxi+ + {\-w) xi, wyi.+ + (1 - и,) у/) -

-a-614+- cji/l+l. (6.4.15)

Аналогичным образом записывается и оценка Еу. Отметим, что в данном случае мы находим отклонение в значениях функций f п g, но не в переменных состояния х п у.

> Обычно норму II II определяют так: !и| = тах\ш\. - Прим. ред.

И^.-. /) = -irH ,-24-f x/,]--i-/t2-g-(, t.), (6.4.16)

тде ge(2, i, 2,+i). Для аппроксимации четвертой производной можно воспользоваться разностной формулой вида

/)-И-2-41-, + б4-4!-ь, + Ы (6.4.17)

TipH 1 = 2,3, п - 2. Тогда оценку погрешности в направлении оси 2 можно записать в виде

£,= max Л-х12-Щ^г + Щ-Щ+г+хи2\- (6-4.18)

i - 2 tl - 2

Каждую из величин £ = тах(£х, £у), приведенных выше, можно использовать для регулировки шага по времени т. С этой щелью обычно выбирается некоторое характерное значение е, причем шаг т уменьшается, если £ > е. Если же £ существенно меньше, чем е, то шаг t увеличивается.

Мы ограничились здесь лишь идейной стороной алгоритма изменения шага т. Очень вероятно, что действительная величина погрешности решения будет заметно отличаться от е. Поэтому если мы хотим удостовериться в том, что заданная точность обеспечивается, следует провести вычисления при двух различных значениях е (отличающихся, например, на порядок) и сравнить полученные результаты. Таким образом, однако, можно оценить лишь влияние погрешностей аппроксимации в направлении оси t . Влияние погрешности аппроксимации в направ-.лении координаты z и адаптационные алгоритмы для выбора шага h будут обсуждаться в следующем пункте.

A.S. Автоматическое изменение шага h

В этом пункте мы исследуем проблему автоматического изменения шага h на равномерной сетке. Для решения вопроса о том, достаточно ли велико п (т. е. достаточно ли мал шаг h), при заданной характерной величине погрешности е, будем искать оценку погрешности аппроксимации в направлении оси 2. В принципе мы могли бы провести вычисления для момента tj дважды - с шагом h и 2h - и результаты сравнить. Такой численный подход не был бы, однако, достаточно простым и эффективным.

Другой, более подходящий путь заключается в том, чтобы использовать выражения для остаточного члена в соответствующих разностных формулах. При замене второй производной в (6.4.3) имеем

Аналогичное выражение получается и для £у (при i = 1 и i = - п-\ для четвертой производной можно использовать асимметричные формулы). Обозначим £ = тах(£х,£у).

Тогда для изменения (регулировки) п можно использовать следующую стратегию:

1. £ > е; шаг h слишком велик, и его следует уменьшить вдвое. Значения решения в узловых точках, которые появятся между исходными узловыми точками, после уменьшения шага, можно определить с помощью интерполяции, например, положить

и./2 = Т^[-1-. + 9А^1 + 94+.-Ы.=1. 2-----п-2,

42 = тИ + 64-4]. (6-4.19)

4-1/2 = Т + 64-1-4-2]-

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

2. £ < е/(о, где значение ш, как правило, выбирается между 4 и 8. Шаг h излишне мал и может быть удвоен. При этом узловые точки с нечетными индексами 1, 3, 5.....п - 1 выбрасываются, оставшиеся точки перенумеровываются, шаг h удваивается, а п уменьшается вдвое.

3. В остальных случаях шаг h остается неизменным и проводятся дальнейшие вычисления.

При формулировке алгоритма удобно задавать максимальный и минимальный допустимый шаг h.

6.4.4. Адаптивная неравномерная сетка

Пусть у нас имеется некоторая, вообще говоря, неравномерная сетка узловых точек 2о = 0, Zu 2 ==1. Тогда вторую производную в направлении z в точке (г,-, tj) можно заменить трехточечной разностной формулой

~(.-.-,)(.-ы-.-1)

-(.-.-1)к.-.) +(...-J(...-.)

При этом разностные уравнения (6.4.3) изменяются только за счет подстановки формулы (6.4.20) вместо простейшей фор-

.(1) -Ci- -) J-ZilfbLZifblr

7-1/2 4(г,-г, ,){г, ,-г^ )1-2- К-.--з)

ВО втором случае-интерполированием по узловым точкам

Zi-i, Zi, Zi+i,

v(2) f + l , i + l i-l

i-U2 Цг,-г, ,) -1+ 4(г,,-г,)

- -.-i) (6.4.22)

4(2: + l-Z,- ,)(z, + ,-Zj)

Для оценки погрешности аппроксимации на интервале {Zi-i, 2/) используется величина

£ = U<-i,/2-xP2,/2l; (6.4.23)

£у определяется аналогично. Тогда £ = max (5, 5J,). Далее сетка узловых точек видоизменяется с помощью следующего алгоритма (здесь снова е - заданная характерная величина погрешности) .

1. £ > е. Тогда между 2; и 2j i мы добавляем еще одну узловую точку (2j i -j- Zi) /2. Приближенное значение решения в этой точке находим посредством интерполяции по четырем соседним узловым точкам 2, 2, Zi-i, Zi, Zi+i. Проделаем эту операцию для i= 1,2, п. Далее перенумеруем узловые точки (добавим новые точки с сохранением порядка), после чего продолжим вычисления на построенной более густой сетке.

2. £< е/(о для i = k л одновременно для i = &-J- 1. В этом случае мы отбрасываем узловую точку Zk и вновь проводим

д^х 1

) g~-r(*i i - 2-*;Формула (6.4.20), очевидно, переходит в простейшую при - г. = г^, - Zj = h. - Прим. ред.

мулы) ДЛЯ -gp-. Все остальное остается неизменным, включая

вычисления нового слоя Погрешность аппроксимации

в направлении оси z оценивается на каждом подынтервале {Zi-u Zi) отдельно, Айгенбергер и Батт [6.29] предложили следующий подход. Приближенное значение решения при / = tj и z = {Zi-i + Zi)/2 отыскивается двумя способами (ниже индекс / опущен). В первом случае мы получаем его интерполированием по узловым точкам 2, 2, 2,-1, Zi,

соответствующие вычисления £ для соседних узловых точек. Затем перенумеровываем узловые точки и продолжаем вычисления. При этом значение ш выбирается обычно в диапазоне 10-20.

3. В остальных случаях сетка не изменяется. При этом бывает удобным задавать максимальное и минимальное расстояние между двумя соседними точками.

Для оценки погрешности аппроксимации на неравномерной сетке можно было бы использовать также соображения, изложенные в предыдущем пункте для случая равномерной сетки. В последнее время появился ряд работ по этой теме, см., например, сборник [6.35].

6.4.5. Метод прямых

С методом прямых мы уже встречались в п. 6.3.3 (см. уравнения (6.3.31)). Его применение оказывается эффективным в тех случаях, когда число узловых точек сравнительно невелико. Последнего можно добиться увеличением порядка аппроксимации. В данном пункте на примере системы типа реакция-диффузия (4.3.7) будет продемонстрировано использование этого метода для случая, когда пространственные производные заменяются разностными отношениями с пятью узловыми точками. При этом порядок аппроксимации по сравнению с трехточечной схемой (6.3.31) возрастает с 0(/t2) до 0(h*). Обозначим д:,(/) x(zi, i), yi{i) та y(Zi, t) и рассмотрим равномерную сетку 2,= = ih, 1 = 0, 1, h=l/n. Уравнение (4.3.7а) заменяется

системой уравнений

124 i-Xi-2+lQXi i - 30Xi + l6Xi+i - Xi+2) + f{Xi, Уд,

i = 2, 3, п-2, (6.4.24a)

4r=lW - 20 + + -xd + f {хг, Уг), (6.4.24b)

+ 4a; 3 - Xn-,) + f (Xnu Уп-i). (6.4.24c)

В случае граничных условий 1-го рода вида (4.3.9), мы имеем Xo{t) х, Xnit) = X. Для уравнения (4.3.7b) соответствующая замена строится совершенно аналогично. Таким образом, мы получаем систему обыкновенных дифференциальных уравнений, которую можно решить каким-либо из методов, описанных в § 5.7, например методом Рунге-Кутты.

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

д^у , а ду дг г дг

J]Aayiri,t). (6.4.25)

1 = 1

Коэффициенты Aif находятся из требования, чтобы левая часть соотношения (6.4.25) точно равнялась правой для функций г/=1, у = г^, у = г*, г/ = г2 -2. Последовательно под-

ставляя эти функции в формулу (6.4.25), получим

2 + f 2г. = 2;Л,Д (6.4.26)

(=1

п

(2п - 2) (2п 3) rf- + (2п - 2) rf- = Arf-.

При каждом фиксированном / соотношения (6.4.26) представляют собой систему п линейных алгебраических уравнений

ограничение, которое определяется требованием численной устойчивости. Если же мы применяем схему интегрирования с автоматическим изменением шага (например, схему Мерсо-на), то будет выбран относительно короткий шаг интегрирования (~ft2), что, как правило, приводит к большим затратам машинного времени.

Использование метода прямых на очень редкой сетке узловых точек (часто неравномерной) приводит к небольшим системам обыкновенных дифференциальных уравнений, поведение которых можно достаточно эффективно анализировать методами, описанными в гл. 5. При этом часто оказывается возможным перенести полученные качественные выводы на исходные параболические уравнения; в частности, динамическое моделирование (численное решение нестационарных уравнений) можно проводить только в тех областях изменения параметров, где в результате анализа с помощью метода прямых следует ожидать появления каких-либо интересных эффектов.

Продемонстрируем указанный подход, построив аппроксимацию решения задачи 16 с помощью метода прямых [6.10]. Зададим граничные условия (Р16-13) при Ни->-схэ, Sh->-oo и выберем п узловых точек гьГг, Гп=1. Предположим теперь, что пространственный дифференциальный оператор в уравнении (Р16-11) заменяется в узловой точке г,- некоторой линейной комбинацией значений решения в узловых точках, а именно

относительно неизвестных коэффициентов Ац. Найдем решения этих систем при / = I, 2, ..., п - I (в случае / = п это оказывается излишним, поскольку в точке г = I у нас задано у = I).

Используя аппроксимацию (6.4.25) для у и Q, мы получаем из исходных уравнений (PI6-I1), (Р16-12) систему обыкновенных дифференциальных уравнений для функций y,{t) = у{г1, t), е/ = в(г„/), /= 1,2 п -I:

= Л .в, + уРФ^г/,ехр^

1 + Q./y 9,

1 = 1

+ e,./Y

(6.4.27) (6.4.28)

(напомним, что yn{t) = 1 и вп() = 0).

Положим теперь п = 2, т. е. выберем лишь одну внутреннюю точку г\, а также точку г? = 1. Система уравнений (6.4.26) для /= 1 приобретает вид

Л.+ .1 = 0. 2(1 +а)

Л„, + г2Л„ = 2(1 + а),

= Р, = -

2(1 +а) г?-1

(6.4.29)

Дифференциальные уравнения (6.4.27), (6.4.28) принимают в этом случае вид

Lw = p(t/,

d&

1)-Ф^,ехр

l-f-Oi/Y

(6.4.30)

= рв, + уРФ^,ехр^

(здесь уже использованы условия y{t)=\, Q{t)Qy

Весьма существенным является в данном случае выбор координаты Гь Вилладсен и Стюарт [6.30] рекомендуют выбирать в качестве ri нуль подходящего ортогонального полинома. Рекомендуемые ими значения (нули неких полиномов Якоби) даются следующей таблицей:

Для стационарного решения системы (6.4.30) имеет место соотношение

e, = Yp(l - i/i), (6.4.31

аналогичное формуле (Р16-15). Подставляя его во второе из уравнений (6.4.30) (при Oi = 0), получаем

р(1-1)-Ф^,ехр^&11 = 0.

(6.4.32)

В табл. 6.12 приведена зависимость г/, от Ф, подсчитанная с помощью метода отображения параметра Ф для последовательности значений ух. Сравнивая ее с табл. 6.9, можно заметить качественное совпадение зависимостей решений от параметра (г/1 соответствует значению г = 0,4472!). При р = 0,05 зависимость однозначна, при р 0,25 у нас возникают кратные решения, а при р = 0,4 в некотором диапазоне изменения параметра Ф существуют три стационарных решения данной задачи.

Таблица 6.12. Зависимость решений уравнения (6.4.32) от параметра Ф для задачи 16, у = 20, а = О, р = -2,5.

Ф

Р=0,05

Р = 0,2

Р = 0,25

Р=0,4

0,95

0,85

0,75

0,65

0,55

0,45

0,35

0,25

0,05

0,354 0,501 0,617 0,716 0,807 0,977 I.I48 1,338 1,573 1,722 1,908 2,153 3,084 4,379

0,329 0,433 0,496 0,538 0,567 0,603 0,626 0,649 0,682 0,707 0,743 0,796 1,032 1,396

0,321 0Д13 0,463 0,491 0,507 0,519 0,520 0,522 0,532 0,545 0,565 0,597 0,756 1,011

0,298 0,359 0,377 0,377 0,368 0,340 0,311 0,288 0,274 0,271 0,272 0,280 0,336 0,439

Динамическое поведение системы (6.4.30) может подсказать нам, как будет вести себя исходная система двух дифференциальных уравнений с частными производными, динамическое моделирование которой требует больших затрат машинного времени. Выбор п = 3 может дать нам более точные результаты. Более подробно результаты исследования задачи 16 приведены в работе [6.19], а для других задач -в работе [6.4].

6.5. ПЕРИОДИЧЕСКИЕ РЕШЕНИЯ В РАСПРЕДЕЛЕННЫХ СИСТЕМАХ

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

6.5.1. Вычисление периодических решений и процесс продолжения

Периодическое (во времени) решение системы типа реакция-диффузия (4.3.7) с заданными граничными условиями по определению удовлетворяет тождествам

xiz,t + T)x{z,t), y{z,t + T)y(z,t) при 2е[0, 1]. Здесь Т - период данного решения.

f=0 0.78 1.86 2.А6 2.82 3,12

(6.5.1)

Рис. 6.13. Периодическое решение задачи 11 при ГУ1; А = 2, В = 5,45, 0 = 0,008, Dy = 0,004, I, = 0,75; период 7 = 3,12.

Простейшим подходом для нахождения устойчивого периодического решения является использование какого-либо из методов динамического моделирования (см. § 6.4) в течение достаточно большого промежутка времени t, за который решение системы стабилизируется, становясь (приближенно) периодическим. Эта стабилизация может зависеть, конечно, от начальных условий. В качестве примера того, как может выглядеть периодическое решение в распределенной системе, на рис. 6.13 при-