Численное интегрирование бигармонического уравнения в квадратной области Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2013. Вып. 1

УДК 519.635.1

В. И. Ряжских, М. И. Слюсарев, М. И. Попов

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ БИГАРМОНИЧЕСКОГО УРАВНЕНИЯ В КВАДРАТНОЙ ОБЛАСТИ*)

1. Введение. Бигармонические уравнения, представляющие собой дифференциальные уравнения в частных производных 4-го порядка, играют важную роль в механике сплошных сред при моделировании поведения упругих пластин под нагрузкой [1], а также в гидродинамике малых чисел Рейнольдса [2].

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

В [6] впервые методом суперпозиции решений Навье и Леви получено аналитическое выражение для искомой функции в виде рядов, однако ценность его невелика, ввиду плохой сходимости этих рядов в области угловых точек из-за наличия в коэффициентах гиперболических функций, существенно лимитирующих число слагаемых при вычислениях. Данное обстоятельство стимулировало исследования по поиску способов улучшения сходимости аналитического решения с использованием асимптотического закона [7], а также по выбору альтернативной системы ортогональных функций [8-10], упрощающих вычислительный процесс.

Параллельно стали разрабатываться эффективные численные схемы интегрирования краевых задач для бигармонического уравнения [11-14], основным недостатком которых является трудность достижения требуемой степени детализации области интегрирования из-за их неявного характера, что предполагает решение систем линейных уравнений очень высокого порядка с неразряженной матрицей. Разрешение такой коллизии может быть осуществлено, например, путем симбиоза маршевых конечно-разностных схем и метода установления.

Ряжских Виктор Иванович — доктор технических наук, профессор, заведующий кафедрой высшей математики Воронежского государственного университета инженерных технологий. Количество опубликованных работ: 82. Научные направления: теоретические основы явлений переноса и вычислительная математика. E-mail: kafvm@vgta.vrn.ru.

Слюсарев Михаил Иванович — кандидат технических наук, доцент кафедры процессов и аппаратов химических и пищевых производств Воронежского государственного университета инженерных технологий. Количество опубликованных работ: 32. Научное направление: моделирование задач тепломассообмена. E-mail: kafvm@vgta.vrn.ru.

Попов Михаил Иванович — аспирант кафедры высшей математики Воронежского государственного университета инженерных технологий. Научный руководитель: доктор технических наук, проф. В. И. Ряжских. Количество опубликованных работ: 3. Научное направление: краевые задачи для уравнений в частных производных (численные методы). E-mail: mihail_semilov@mail.ru.

*) Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 10-08-00120-а).

© В. И. Ряжских, М. И. Слюсарев, М. И. Попов, 2013

2. Постановка задачи и построение численной схемы. В связи с этим рассмотрим краевую задачу на квадратной области С = [0,1] х [0,1] для неоднородного бигармонического уравнения

д4Ф(х,у) , ^д4Ф(х,у) , д4Ф(х,у) _ ^

дх4 дх^ду! + ду4 - Р[Х'У)> [1)

Ф(0, у) = Ф(1, у) = Ф(х, 0) = Ф(х, 1) = 0, (2)

ЭФ(0,2/) = дФ(1,у) = дФ(х,0) = дФ(х, 1) = 0

дх дх ду ду '

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

Применим процедуру установления, считая, что решение системы (1)-(3) можно получить из решения системы

дФ(х,у,т) д4Ф(х,у,т) д 4Ф(х,у,т) д4Ф(х,у,т)

+—кг1+"^Ш2+—^ = &

Ф(х,у, 0) = 0, (5)

Ф(0, у, т) = Ф(1, у, т) = Ф(х, 0, т) = Ф(х, 1, т) = 0, (6)

дЩ0,У,т) = 0Ф(1 ,у,т) = дЩх,0,т) = 0Ф(х,1,т) = 0

дх дх ду ду '

когда т ^ ж, здесь т - фиктивная переменная (аналог времени) [15].

Разобьем область решения [0,1] х [0,1] (в общем случае - произвольную прямоугольную область) системы (1)-(3) равномерной сеткой с шагами Ах и Ду, 1дАх<Ау = {х, у-) = (гДх,] Ду), г = 0,.., п, ] = 0,.., т}. Вместо функции непрерывного аргумента на С будем рассматривать функцию дискретного аргумента Ф(а^, Уу)&х,Ау Значения функции в узле будем обозначать Ф^-. Выберем разбиение таким образом, чтобы в него попадали точки границы области.

Конечно-разностную аппроксимацию производной дФ(х,у,т)/дт в произвольном узле (г,з) находим применением конечно-разностного оператора 1-го порядка

дФ(х у т) Фк+1 - Фк ■

ТФ = ^ = —- + °(Л

дт Л т

Для построения конечно-разностного аналога частных производных д4Ф(х,у)/дх4 и д4Ф(х,у)/ду4 используем центрально-разностный оператор 2-го порядка, примененный дважды по соответствующей переменной [16]:

д4Щх,у,т) Ф^2,3 - 4Ф?+и + 6Ф^ - 4Ф,-и + В^ - -- -- + °((Л )>

г)4Ф(т у т) Фк •, о - 4Фк •, 1 + 6Фк • - 4Фк ■ 1 +Фк • о ВзФ = ь= '3+ (д ;)4 + + 0((Д уП

Аппроксимацию смешанной производной определим с помощью пятиточечного приближения второй производной [16] по х и по у

- 16*?+2,— + *?+2,—) - 16(Фк+1,^+2 - 16фк+1,д"+1 + 30*?+1- 16*?+!,— + + *?+1,—) + 30(*?,— - 16*?— + 30*,. - 16*?,— + *?,—) -

- 16(*?-1,— - 16*?-1— + 30*?_1,.. - 16*?_1,— + *?— 1—2) +

+ (*?— 2,— - 16*?— 2,— + 30*?—2,. - 16*?—2,. —2 + *?—2,. —2)) + 0((Д х)2 (Д у)2).

Обозначим В* = В1* + В2* + В3*. Заменяя в уравнении (1) частные производные конечно-разностными аналогами, а известную функцию ее сужением на узлы сетки и отбрасывая слагаемые более высокого порядка малости, получим уравнение

Л г (Л ж)4

2

+ 144(Лх)2(Л^((Ф^+2 " 16Ф»Ч-2,д-+1 + 30Ф-, - 16*?+2— +

+ *?+2,. —2) - 16(*?+1,.+2 - 16*?+1— + 30*?+1,.. - 16*?+1,— + + *?+1,—) + 30(*?,— -16*?— + 30*?,.. - 16*?,..—1 + *?—) -

- 16(*?— 1,— - 16*?—1,— + 30*?—1,.. - 16*?—1,..—1 + *?—1,..—2) + + (*?— 2,— - 16*?— 2,— + 30*?—2,. - 16*?—2,. —2 + *?—2,.. —2)) +

+-^-= (8)

которое можно переписать следующим образом:

, * ? „ . _ 4* ? _ . + 6* ? ■ - 4* ? . + *? о .

V и V (Дж)4 +

2

+ 144(Лх)2(Л^((Ф^+2 " 16Ф"г + + 1 + " 16Ф*+2— +

+ * ?+2,. — 2) - 16(*?+1,— - 16*?+1,— + 30*?+1,. - 16*?+1,. —1 + *?+1,— ) + + 30(*?,.+2 - 16*?— + 30*?,.. - 16*?,.. —1 + *?. — 2) -

- 16(*?—1,.+2 - 16*?—1,— + 30*?—1,. - 16*?—1,. —1 + *?—1,. —2) + + (*?— 2,— - 16*?— 2,— + 30*?—2,. - 16*?— 2,. — 2 + *?—2,. —2)) +

, Ф?,д-+2-4Ф?,— + 6Ф?|Д.- 4Ф?—+Ф?|Д._2 Ч

+ (Л у)4 + ^

Начальное условие (5) на сетке приобретет вид

Ф°—0, г,] = 0Я (10)

Дискретный аналог граничных условий, накладываемых на саму функцию, получим сужением условий (6) на граничные узлы сетки

Фо,— Фп,— Ф?,о = Ф?,п = 0, к = 1,2,..., г = 0/п, ] = 0~п. (И)

Аналог граничных условий, накладываемых на частные производные 1-го порядка (3), находим из соотношения для конечно-разностного аналога первой производной

фк - фк . - фк . - -

1,3 0,_7 __п—1 ^ __г,0 _ г,п г,п— 1

= 0, г,з = 1,п — 1.

Л х Л х Л у А у

Учитывая (10), имеем

фу = Фп-1,3 = 41 = Чп-1 = * = 1. 2, * = Т^Х = Т^Т (12)

Таким образом, непрерывная задача (4)-(7) заменена конечно-разностной схемой (9)-(12).

На множестве сеточных функций {и(х^,уз)дх,ду} введем Гильбертово пространство Н дх,Ау, которое для каждого разбиения представляет собой вещественное пространство (п + 1) х (т + 1)-мерных векторов. Чтобы учесть граничные условия (11), (12), потребуем их выполнение для каждой функции из Ндх,ду. Зададим скалярное произведение на Ндх,ду в виде [17]

п — 2 т — 2

(и, у) = Л х Л у.

г=2 г=2

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

3. Вычисление погрешности аппроксимации бигармонического уравнения разностной схемой. Рассмотрим уравнение (8) в операторной форме

Ти + Ви = —/. (13)

Для оценки его точности образуем разность г = и — у, где и - решение задачи (9)-(12), а V - решение задачи (4)-(7) [17]. Подставляя и = г + у в (13), получим для г задачу

Ти + Ви = Т (г + у)+ В(г + у) = (Т + В) г +(Т + В)у = —/, Вг = —ф.

Во внутренних узлах сетки и для г выполнены граничные условия (10)—(12), где ф = (Т + В)у + / - погрешность аппроксимации задачи (4)-(7) схемой (9)-(12). Так как + Д2-у + / = 0, то

дV ду

ф = (Т + В)у + /- — + А\-/=(Т + В)у-(— + А2У).

Вычислим погрешность аппроксимации почленно, используя разложение по формуле Тейлора с остаточным членом в форме Лагранжа, в окрестности г,з, к-го узла:

„ , N д4у(х,у) Л х2 д6у, . ,

ВМх, у, т)--= __[д х{% + в), % А у, к А г], 0 < в < 1,

д4у(х, у) Л у2 д6у

В3у(х,у,т)--= -¿-д^Ь А х,Ау{3+в),кА г], 0 < в < 1,

„ , ч д4у(х, у, т) 17 (Л х6 д8у , ,

ВМ-, У, г) - 2 ¿2ду2 = — ^ [Д + 913 л у, к А г] +

л у6 д8у , , \ 68 (Л х5 д8у ^

+ + Т]) " 189 + Л УЬ +

Д у5 д8« . , \ 167 (

167 ( 4 д8г

Д х дхду7

дх6ду2

д 8у ч

х [Д х(г + 6»), Д у+ а), к А т]+ Д у4^^ 6 [Д х(г + 6»), Д + а), к А т]J -

дх2ду6

68 /

д 8г

- ^(Д ^ Ау^-^[Ах(г + в),А у{3+а),к А т] +

д8у \

+ Ау3 [Д + +

85 д 8«

+ — Д х2 А У2 дх4ду4[А х{г + в),А у{з + а),к А т]; 0 < 0 < 1, 0 < а < 1, „ д«(х, у,т) Д т

д2

V г ,

Тч)--э = 2 э 2[гАж,гДу,Дг(г + р)], 0 < р < 1.

На равномерной сетке Д х =Д у = Н. Обозначим

д6«(х, у, т)

Мх = тах

Мз = тш

(х,у,т)еах[ 0,4]

М4 = тах

(х,у,т)еОх[ 0,4]

дх6 д8у(х, у, т)

, М2 = тах

(ж,у,г)еСх[0,4]

д6«(х, у, т)

ду6

дхтдуп д8«(х, у, т)

дхтдуп

М = тах

(х,у,т)еОх[ 0,4]

т, п = 1, 3, 5, 7, т + п = 8,

, т,п = 0, 2, 4, 6, 8, т + п = 8, 1 д2у(х, у, т)

2 дт2

где -^время установления. В силу симметрии оператора М1 = М2, обозначим М1 М1/3, М2 = (5398/945)М4 - (1088/189)М3. В результате получим оценку

\ф\ < М Д т + М1Н2 + М2Н4.

(14)

4. Сходимость и устойчивость. Рассмотрим уравнение (8) в операторной форме (13) с начальным и граничными условиями (10)-(12) на равномерной сетке Д х =Д у = Н и Д т = т. Операторы Т : Н^ ^ Н^, В : Н^ ^ НПредположим существование решения задачи (4)-(7), удовлетворяющего уравнению (13). Схема (13) является простейшей двухслойной итерационной схемой. Перепишем ее в канонической форме [17]

—^ + вкик = к = 0,1,..., задан и0,

(15)

здесь А? = Е, В? = В для всех к, ио = 0.

Итерационная схема (15) при любых т точно аппроксимирует уравнение (1) на его решении V. Поэтому разность г? = и? - V удовлетворяет однородному уравнению [17]

£?+1 - г?

+ Вг? = 0, к = 0,1,..., задан г0 = и0 - V.

Решая уравнение (15) относительно г¡?+1, получим

г? +1 = Ягк, Б = Е - тВ,

здесь Б - оператор перехода.

Рассмотрим оператор В. С учетом равномерности сетки

х

1 ( 1 2 17 2 1

/Й \72 ~ 9 + 1 + 12 Щ+2'3 ~ д "¿+2,5-1 + ^ "¿+2,^-2 -

2 32 32 32 2 17

- д "¿+1,5+2 + у "¿+1,5 + 1 - у "¿+1,5 + У "¿+1,5-1 - д "¿+1,5-2 + "¿,5 +2 -

32 49 32 17 2 32

~ У "¿,¿+1 + у "¿,5 - у "¿,¿-1 + "¿,5-2 - д "¿-1,5 + 2 + у "¿-1,5 + 1 -

32 32 2 1 2

- у "¿-1,5 + у "¿-1,5-1 ~ д "¿-1,5-2 + ^ "¿-2,5 + 2 ~ д "¿-2,5+1 +

17 2 1

17 2 1

+ ^ "¿-2,5 ~ д "¿-2,5-1 + ^ "¿-2,5-2^ •

Обозначим В = Ь4В оператор, не зависящий от Ь. Утверждение 1. Оператор В самосопряженный.

Для д о к а з а т е л ь с т в а утверждения необходимо заметить, что матрица оператора В вещественна и симметрична [18].

Утверждение 2. Оператор В положительный.

Доказательство. По определению положительности оператора (Ви, и) > 0 для любых и € Н^, и = 0. То есть

г —2 т — 2 п—2 т—2

г=2 г=2 г=2 г=2

53 53 ВщззЛ х Л у = Ь2 53 53 Ви> з> 0. (17)

Поскольку Ь? > 0, достаточно показать, что ^п=2 Вщ^> 0. Каждое слагае-

мое этой суммы имеет вид

( 1 2 17 2 1 2

"¿+2,5+2 - - "¿+2,5 + 1 + — "¿+2,5 - - «¿+2,5-1 + ^ "¿+2,5-2 ~ д "¿+1,5+2 +

32 32 32 2 17 32

+ у "¿+1,5 + 1 - у "¿+1,5 + "д" "¿+1,5-1 - д "¿+1,5-2 + "¿.¿+2 ~ у "¿,5+1 +

49 32 17 2 32 32

+ у "¿,5 ~ у "¿,5-1 + ^ "¿,5-2 - д "¿-1,5+2 + у «¿-1,5+1 ~ у "¿-1,5 +

32 2 1 2

+ у "¿-1,5-1 - д «¿-1,^-2 + — "¿-2,5 + 2 - д "¿-2,5 + 1 +

17 2 1

17 2 1

+ — "¿-2,5 - д "¿—2,5 — 1 + —"¿-2,5-2^ "¿,5 • (1»)

В силу симметрии выражения (18) каждое произведение координат будет входить в сумму в виде

49 2 49 2

—и^ + 2 + — и1р, (19)

где К - один из 5 коэффициентов Щ-, Покажем положительность выра-

жения (19) для множителей, отличных от 0:

49 2 49 2

У"г,5 + 2Киг,Зи1,Р + ~2~и1,Р ^ Щиг,Зи1,р\ + >

> 49\щзщ,р\ — 2\Кщ,зщ,р\ > 18\щзщ,р\ > 0. Если один из сомножителей равен нулю, то положительность (19) очевидна. Так как и = 0, то хотя бы одно из выражений (19) положительно. Поскольку сумма (17) состоит хотя бы из одного положительного и остальных неотрицательных слагаемых, она положительна. То есть оператор В положительный.

Оператор Б = Е - тВ - линейная комбинация самосопряженных операторов, поэтому он самосопряженный [18]. Из (16) находим \\г?\\. Оценим \\г?\\

1Ы1 < \Т?\\\\го\\ = \\Б\\? \\го\\ , (20)

где Тк = Б? - разрешающий оператор, причем ||Б?|| = \\Б\\?, поскольку Б = Б*. Итерационный процесс сходится, если \\и? - «\\ = \\г?\\ ^ 0, к [3], где \\ • \\ - некоторая норма в пространстве Н^. Поэтому нужно, чтобы \\Б\\ ^ 1. Количество шагов к, необходимых для обеспечения точности е, т. е. \\г?\\ ^ £ \\го\\, посчитаем по формуле

в которой величина ln(1 / У Б ||) - скорость сходимости. Чем меньше ||Б||, тем меньше n(e), соответственно скорость сходимости больше. Поскольку ||Б|| зависит от т, выберем т из условия минимума ||Б||. Так как Б - самосопряженный оператор, то

ЦБЦ = max\pk\,

k

где yUfc - собственные значения. Из определения Б (16) следует, что ри = 1 — jjzХк, ^к -собственные значения оператора B. В результате получаем задачу минимакса

min max 11 — TAk/НА\.

т>0 k

Лемма [17]. Если 0 < Amin ^ А ^ Amax, т > 0; то

I л I Amax Amin 2

mm max |1 — тА| = ----= po при т = tq =

i i Л..1. А. Л. Л. A. XV.A'd I J- JS 0 I иI >0

r>0 Ae[Amin,Amax] Amax + Amln AmaX + A

mm

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

A — A . 2h4

II «и Amax Amin /oo\

II II = T-—T-=Po при T = TO = --—-. (22)

Amax + Amln Amax + Amln

Таким образом, получено не только оптимальное значение У Б||, но и выражение для оптимального значения шага (22). Поскольку вычисление собственных значений оператора В довольно трудоемкая операция (для сетки 25 х 25 необходимо посчитать собственные значения матрицы размера 441 х 441), используется другой способ вычисления по формулам [19]

Amax = sup (В, u), Amin = inf (В, u)

||и||=1 ll«ll = l

с помощью алгоритмов оптимизации Maple. В табл. 1 приведены расчеты ||Б||, оптимального шага то и количества шагов, вычисленного по формуле (21).

Поскольку ||Б|| < 1 для любых размерностей сетки, схема (15) сходится. Действительно,

u - «|| = Ы < ЦБЦк ||zo| ^ 0, к ^ю.

B = sup Bu

C 1М1с = 1 C

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

\Ы\с = \ич \,

800

а минимальное значение стремится к 0.

Рассмотрим схему (15) как обычную двухслойную схему. Докажем ее устойчивость. Для этого воспользуемся следующей теоремой [17].

Теорема. Для устойчивости схемы (15) достаточно, чтобы выполнялось условие

\\Tnj\\ ^ М1 при любых 0 ^ ] ^ п.

При этом для решения задачи (15) верна априорная оценка

\К+1\\(1) < М1 ^\Ы\(1) + ^ Г \\A-r1 ¡3\\(1) | .

Разрешающий оператор Тп,з = 8п_1Бп_2...8з.

Поскольку оператор перехода схемы (15) постоянный с нормой, меньшей единицы, очевидно неравенство

\\Тп,з\\ ^ \\$\\ = р0 при любых 0 < ] < п. Таким образом, схема (15) устойчива и выполняется оценка

IK+1II < ро ( IM

Итак, с одной стороны, рассматривая схему (15) как явную итерационную схему с постоянным параметром т, мы установили ее сходимость к решению задачи (4)—(7), а следовательно, задачи (1)—(3) и вывели формулы для нахождения оптимального значения параметра т и необходимого количества итераций K. C другой стороны, рассматривая (15) как двухслойную схему в общем виде, мы выявили ее устойчивость и получили оценку для максимального по модулю значения решения. Стоит заметить, что схема (15) однозначно разрешима для достаточно малых h, поэтому ее сходимость следует из устойчивости и выражения (14) для погрешности аппроксимации [20]. Действительно, с учетом (22) погрешность аппроксимации схемы (15), определяемая формулой (14), ф = O(h2).

5. Анализ решений. В виртуальной среде Maple 13 проведены вычислительные эксперименты. Решена задача (1)—(3) с постоянной правой частью, равной —1. В качестве начального приближения взята функция ио = 0. Эксперимент осуществлялся для различных значений размерности сетки и параметров т - мелкости разбиения, е -точности и K - количества итераций. Результаты эксперимента, а также вычисления входных данных представлены в табл. 1, 2 и на рис. 1, 2.

Данные расчетов входных параметров приведены в табл. 1, в которой ^^|| - норма оператора перехода, £ = К2 - точность, К - минимальное число итераций, необходимое для обеспечения заданной точности, то - величина итерационного шага, \итах\ - верхняя априорная оценка максимального по модулю значения решения. Точность была выбрана таким образом, чтобы она совпадала с погрешностью аппроксимации [17].

Таблица 1. Расчет входных параметров

Размерность сетки Мелкость разбиения 11^11 £ К то | ^тах |

11 X 11 0.1 0.995 0.01 895 2.45 • 10~ь 0.0026

15 х 15 1/14 0.9989 1/196 4831 6.09 • Ю-7 0.0029

21 х 21 0.05 0.9998 0.0025 25970 1.43 • 10-у 0.0037

25 х 25 1/24 0.9999 1/576 59762 6.86 • Ю-8 0.0041

Таблица 2. Результаты вычислительных экспериментов

Размерность сетки то К | ^тах | Сходимость

11 X 11 2.45 • 10-ь 895 7.512 • Ю-4 Сходится

11 X 11 2.45 • 10-ь 1000 7.545 • Ю-4 Сходится

11 X 11 2.47- 10~ь 895 0.017 Расходится

11 х 11 1•10~ь 895 6.374 • Ю-4 Сходится

11 х 11 2.45 • 10-ь 500 6.985 • Ю-4 Сходится

15 х 15 6.09 • 10-у 4831 8.97- Ю-4 Сходится

15 х 15 6.09 • 10-у 5000 8.98 • Ю-4 Сходится

15 х 15 6.1 • ю-7 4831 0.016 Расходится

15 х 15 1 • ю-у 4831 5 • Ю-4 Сходится

21 х 21 1.43 • 10-у 25 970 1 • 10~'6 Сходится

21 х 21 1.44 • Ю-7 2000 0.9 Расходится

25 х 25 6.86 • Ю-8 59 762 1.05 • Ю-3 Сходится

О

-0.0002

-0.0004 и(х,у)

-0.0006

-0.0008

-0.0010

В табл. 2 представлены результаты численных экспериментов. В колонке «Сходимость» указано, сходится ли алгоритм к решению задачи. Сравнивая табл. 1 и 2, видно,

и

О

-0.00010 -

-0.00020 -

-о.ооозо -

-0.00040 -

-0.00050 -

-0.00060 -

-0.00070 -

-0.00080 -

-0.00090 -

-0.00100 -

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Рис. 2. Эволюция решения по размерности сети Размерность сетки: 1 - 11 х 11, 2 - 15 X 15, 3 - 21 х 21, 4 - 25 X 25.

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

На рис. 1 изображен график решения u(x, y) для размерности сетки 25 х 25 при оптимальных значениях то и K. Легко заметить, что он симметричен относительно центральной точки x = 0.5, y = 0.5, что является следствием самосопряженности оператора перехода S. В этой точке достигается максимальное по модулю значение функции \umax\. Потому строим сетки так, чтобы центральная точка области попадала в узел разбиения.

На рис. 2 изображена эволюция решения по количеству узлов сетки в плоскости x = 0.5, начиная от 11 х 11 и заканчивая 25 х 25. Хорошо видно, что при увеличении размерности сетки значения решения u(x,y) в тех же узлах все меньше отличаются друг от друга. Это говорит о том, что данный процесс сходящийся.

6. Заключение. Отметим плюсы предложенной конечно-разностной схемы. Во-первых, ее преимущество в том, что она явная, и не приходится решать систему уравнений с неразряженной матрицей, что весьма облегчает вычислительный процесс. Во-вторых, можно заранее задавать точность и в соответствии с ней вычислять оптимальный шаг и количество шагов, необходимых для сходимости алгоритма. В-третьих, явный характер схемы в отличие от неявного позволяет хранить в памяти ЭВМ только один шаг итерации, что существенно экономит память ЭВМ и время расчета. Решение задачи согласуется с известными результатами С. П. Тимошенко [1].

Литература

1. Тимошенко С. П., Войновский-Кригер С. Пластины и оболочки. М.: Наука, 1966. 636 с.

2. Хаппель Дж., Бреннер Г. Гидродинамика при малых числах Рейнольдса / пер. с англ. В. С. Бермана, Г. В. Маркова; под ред. Ю. А. Буевича. М.: Мир, 1976. 630 с. (Happel J., Brenner H. Low Reynolds numbers hydrodymics.)

3. Постнов В. А., Ростовцев Д. М., Суслов В. П., Кочанов Ю. П. Строительная механика корабля и теория упругости: в 2 т. Л.: Судостроение, 1987. Т. 2. 416 с.

4. Слюсарев М. И., Чертов Е. Ю., Ряжских В. И. Аналитическое решение первой тестовой задачи свободной конвекции для кондуктивно-ламинарного режима // Вестн. Воронеж. гос. техн. ун-та. 2010. Т. 6, № 7. С. 165-167.

5. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. 742 с.

6. Бубнов И. Г. Напряжения в обшивке судов от давления воды // Морской сборник. 1902. Т. 312, № 10. С. 119-138.

7. Чехов В. Н., Пан А. В. Об улучшении сходимости рядов для бигармонической задачи в прямоугольнике // Динамические системы. 2008. Вып. 25. С. 135-144.

8. Суслов В. П., Качанов Ю. П., Спихаренко В. Н. Строительная механика корабля на основе теории упругости. Л.: Судостроение, 1972. 720 с.

9. Selvadurai A. P. Partial differention equations in mechanics 2. New York: Springer, 2004. 698 p.

10. Слюсарев М. И., Чертов Е. Ю., Ряжских В. И., Богер А. А. Кондуктивно-ламинарная естественная конвекция ньютоновской тепловыделяющей жидкости в квадратной каверне с постоянной температурой стенок // Вестн. Воронеж. гос. ун-та. Сер. Физика. Математика. 2011. № 1. С. 214-218.

11. Кобельков Г. М. О сведении краевой задачи для бигармонического уравнения к задаче типа Стокса // Докл. АН СССР. 1985. Т. 283, № 1. С. 539-542.

12. Флетчер К. Вычислительные методы в динамике жидкостей: в 2 т. / пер. с англ. А. И. Державиной; под ред. В. И. Шидловского. М.: Мир, 1991. Т. 1. 504 с. (Fletcher C. A. J. Computational techniquies for fluid dynamics.)

13. Алгазин С. Д. Численные алгоритмы классической математической физики. М.: Диалог-МИФИ, 2010. 240 с.

14. Galanin M. F., Milyutin D. S., Savenkov E. B. Development, research and application of a finite superelements method for solution biharmonical equation: preprint Inst. Appl. Math., the Russian Academy of science. М., 2005. P. 1-24.

15. Марчук Г. И. Методы вычислительной математики. М.: Наука, 1980. 535 с.

16. Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная механика и теплообмен: в 2 т. / пер. с англ. С. В. Сенина, Е. Ю. Шальмана; под ред. Г. Л. Подвидза. М.: Мир, 1990. Т. 1. 384 с. (Anderson D. A., Tannehil J. C., Pletcher R. H. Computational fluid mechanics and heat transfer.)

17. Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 552 с.

18. Треногин А. А. Функциональный анализ. М.: Физматлит, 2002. 488 с.

19. Ахиезер Н. И., Глазман И. М. Теория линейных операторов в гильбертовом пространстве. М.: Наука, 1966. 544 с.

20. Филиппов А. Ф., Рябенький В. С. Об устойчивости разностных уравнений. М.: Гос. изд-во теор.-техн. лит., 1956. 171 с.

Статья рекомендована к печати проф. А. П. Жабко. Статья принята к печати 25 октября 2012 г.