CC BY
171
УДК 539.3, 518.6
В.А. Крысько, Л.Ф. Вахлаева, Т.В. Молоденкова РАЗНОСТНЫЕ СХЕМЫ ВТОРОГО И ЧЕТВЕРТОГО ПОРЯДКОВ ТОЧНОСТИ ДЛЯ УРАВНЕНИЙ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
При решении начально-краевых задач для уравнений математической физики разностным методом важным вопросом является выбор порядка аппроксимации по пространственным координатам, а также поиск экономичного алгоритма решения соответствующих систем разностных уравнений. В работе построены схемы повышенной точности для многомерного слабонелинейного уравнения теплопроводности, волнового уравнения, уравнения колебания пластинки, уравнения движения пластинки в кинематической гипотезе Кирхгофа. Проведено сравнение со схемами второго порядка точности на модельных задачах. Выявлен экономичный алгоритм решения нелинейных разностных уравнений для каждой задачи.
V.A. Krysko, L.F. Vakhlaeva, T.V. Molodenkova DIFFERENCE SCHEMES OF SECOND AND FOURTH ORDERS OF ACCUARANCY FOR EQUATIONS OF MATHEMATICAL PHYSICS
At the decision of initial value or boundary value problems for the equations of mathematical physics by difference method the important question is the choice of approximation order on spatial coordinates, and also search of economic algorithm of the decision of corresponding systems of difference equations. Schemes of the raised accuracy for many-dimensional weakly nonlinear equations of heat conductivity, the wave equation, the equation of fluctuation of a plate, the equation of movement of a plate in Kirhgof's kinematics hypothesis are constructed in this work. Comparison with circuits of the second order of accuracy on modelling problems is carried out. The economic algorithm of decision of nonlinear difference equations for each problem is revealed.
При численном решении начально-краевых задач для уравнений математической физики разностным методом возникает проблема выбора того или иного порядка аппроксимации разностной схемы. Чем выше порядок аппроксимации по пространственным координатам, тем меньше порядок системы разностных уравнений. При решении нестационарных задач аппроксимация второй производной по времени в уравнении и первой производной в начальном условии имеет только второй порядок. Поэтому для получения требуемой точности приходится выбирать достаточно мелкий шаг по времени. Повысить порядок аппроксимации по времени возможно, если применить метод прямых в сочетании с разностным методом по пространственным координатам, а для решения системы дифференциально-разностных уравнений использовать метод Рунге-Кутта четвертого порядка. В этом случае начальные условия (dwldt=ty2(x,y)) аппроксимируются точно. Исследования разностных схем проводи-
лись на модельных задачах для следующих уравнении: одномерного и двумерного волновых уравнений, уравнения колебания пластинки, уравнения движения пластинки в кинематической гипотезе Кирхгофа и слабонелинейного многомерного уравнения теплопроводности. По каждой задаче был найден оптимальный алгоритм.
1. Начально-краевая задача для одномерного волнового уравнения: д Т д Т
ди = ди + /(X, 0,0 < X < 1,0 < і < Т, (1)
д і д х
и(х,0) = ио (х), 50(Х,0) = и1 (х), 0 < х < 1 , (2)
ді
и(0,і) = д1(і), и(1,і) = д2(і), 0 < і < Т . (3)
Введем сетку тоИт = (х; = іИ, і = 0,Ы, И = 1/N; іі = jт, ] = 0,М, т = Т/М} и сеточную функцию у/ = у (хі ,і}.) . Задаче (1)-(3) сопоставим разностную схему с весами
Уіі =ЛУ(а) +Ф, (х,і) є®ит , (4)
У0 = U0(^ y'= У0 + 2(U0(хг)+f(x,,0))+0(т2), i =1N-1 ,
У00 = дl(ti), yN = Д2(tj), j = 0M , (5)
где
у-- = (У/+1 - 2yi + У/-1 Vт2 , ЛУг = Ухх = (Уі-1 - 2Уг + Уг+l Vh2, y(a) = ay1+1 + (1 - 2a)y1 + ay1-1, a = const > 0.
Вводя функцию погрешности решения z=y-U, подставляя в (4) вместо y=z+U, получим разностное уравнение относительно погрешности z: z-tt = Лz(a) + у, (x,t) єшйт, где
у = ф + Ли(a) - U-tt является погрешностью аппроксимации разностной схемы (4) на решении
U(x,t) исходной задачи (1). Если схема аппроксимирует исходную задачу и устойчива, то она сходится, и порядок точности схемы совпадает с порядком аппроксимации [1]. При a=0 схема (4) - явная, она условно устойчива при T<h, погрешность аппроксимации y=0(T2+h2), если фі = f (хг ,tj). При a>l/4 схема неявная, безусловно устойчива, погрешность аппроксима-
2 2 1 h2 i r/ \ h2 д2f
ции y=0(h +т ), если ф1 = f (хг ,ti). При a = a, =------------- и (f)J = f (Хг, t i ) +-----—
V 7 V, V . 4(1 -Є) 12т^ 12 дх2
погрешность аппроксимации y=0(h4+T2), это схема повышенной точности. Для определения yj+1 получаем из (4)-(5) трехточечную разностную краевую задачу
aY2 (yd1 + У?!) - (1 + 2aY2)У/+1 = -F , г = 1N -1, У0+1 = дГ1, УЇЇ1 = д2+1 ,
т h ’
F, = (2 у, - y/-1 + т 2(1 - 2a) Л у1 +aт2Л у1-1 +т2ф, у = т
которая решается методом прогонки.
Исследования на модельных задачах показали, что оптимальным алгоритмом среди неявных схем является схема повышенной точности 0(И4+т2), т.к. шаг по пространственной координате в 10 раз больше, чем для схем 0(И +т ), что приводит к уменьшению в 10 раз порядка системы разностных уравнений, а, следовательно, к уменьшению затрат машинного времени.
Модельные задачи соответствуют точным решениям:
1. и ( х, і) = х2 і2 + 5 х +і;
2. и (х, і) = (бій пх)і + х2;
3. и (х, і) =х4і + хі3 + 2;
в области В = {0 < х < 1,0 < і < 10} .
(к)
Результаты исследований приведены в табл. 1, где 1ЫГ' = тах у/- и/
И,т
погреш-
ность решения на сетке ш*т, £=1,2,3 - номер задачи.
Таблица 1
Влияние порядка аппроксимации на шаги сетки Л и т
Порядок аппроксимации 0(Л4+т2) 0(Л2+т2)
с с* 0.5 1/3
Л; т 0.1 ; 0.01 0.01; 0.01 0.01; 0.01
N1 И,Т 0.037 0.040 0.039
N1ИТ 0.025 0.031 0.030
N1И3Т 0.042 0.045 0.043
Лі; ті 0.05; 0.001 0.001; 0.001 0.001; 0.001
ИИ1,т1 0.009 0.012 0.012
Применим теперь к задаче (1) метод прямых
^ = Р, 1 < I < N -1, йг
(6)
йР
~Г = Ухх,г +Ф(X,г), X . аг
Порядок системы дифференциально-разностных уравнений равен 2^-1). Начальные условия (2) задаются точно.
У (0) = и0( х), Р (0) = ^( X), 1 < г < N -1. (7)
Краевые условия:
= Д:(г;), =^2(г;). (8)
Решая систему дифференциально-разностных уравнений (6)-(8) методом Рунге-Кутта 4-го порядка, получим окончательно погрешность аппроксимации у=0(Л4+т2), если Ф=(х, г);
у=0(И4+т2), если ф = / (х, і) +
И2 д 7
12 дх2 '
Решение системы дифференциально-разностных уравнений (6) находим по формулам Рунге-Кутта 4-го порядка
УҐ = У/ +т( + 2к2 + 2кз + к4), і = 1,N -1 , (9)
6
где У1 = (у/, у2,..., у]ы-1,Р11 ,Р/,..., Р^) - искомый вектор решения системы (6), записанной в векторном виде:
йУ
= Р (г у) . (10)
аг
В (9) кг вычисляются по формулам
£3 = Р| г1 + ^,у1 + -Гк2 | , £4 = Р(г1 +т,у1 +т£3).
Методы Рунге-Кутта являются явными (о=0) и одношаговыми, счет идет по явным формулам. На практике рекомендуется проводить расчеты на нескольких сгущающихся сетках. Если при сгущении сетки решение мало меняется, то требуемая точность достигнута.
Для решения модельных задач 1, 2, 3 применим метод Рунге-Кутта (9) с аппроксимацией по пространственной координате 0(к4) и 0(к2). Результаты приведем в табл. 2. Сравнение табл. 1 и 2 показало, что наиболее экономичным алгоритмом является метод прямых с погрешностью аппроксимации у=0(т4+к4), т.е. метод Рунге-Кутта 4-го порядка 0(т4) с аппроксимацией по х 0(к4).
Таблица 2
Влияние порядка аппроксимации в методе Рунге-Кутта на число уравнений
Порядок аппроксимации Рунге-Кутта 0(Л4+т4) Рунге-Кутта 0(Л2+т4)
Л; т 0.1; 0.001 0.01; 0.001
Число уравнений 2(^-1) в (10) 18 198
N11д 0.009 0.0097
1И1(2) 0.009 0.009
N1 к! 0.009 0.011
2. Начально-краевая задача для двумерного волнового уравнения:
д 2и и
ди- = ьи + /(х,г), Ъ = {0 < х = (х1, х2) < 1,0 < г < Т} = Ъ + Г , (12)
дг
и(х,0) = и0(х), ди(х,0) = и1(х), 0 < х = (х1,х2) < 1,
дг (13)
и|Г = |д(х,г), 0 < г < Т,
д2и
где ьи = ьи + ь2и, ьаи = ^, а = 1,2.
д ха
а
Введем сетку токт = {х, = (х1(г),х21)); г,, = 0,N, к = 1/N; гк = кт, к = 0,М, т = Т/М}, и в узлах сетки сеточную функцию у£ = у (гк, ,'к, гк ) .
Поставим в соответствие задаче (12)-(13) разностную схему с весами [1]
У-П = (Л1 +Л2)У(a) +Ф, (х-) Є «йт, х = (xl(г),х21 ^ тойт = «йт +Yйт , (14)
у(х,0) = U0(х), У-(х,0) = U(х) + 0,5т(LU(х,0) + f (х,0)) + 0(т2),
У 1„ =д(х-). (15)
где ЛaУ = y*a*a, a= 1,2, ф= f (хЛ), У(a) = ayk+1 +(l - 2a)yk + ayk-l.
По аналогии с одномерным волновым уравнением вычислим погрешность аппроксимации
у = (Л1 + Л 2 )U (a) - U-t + ф .
При a=0 получим явную схему, она условно устойчива при т< h/V—, у=0(й2+т2),
f 2 YV—
z=0(h2+т2) [1]. При a<l/4 схема неявная, условно устойчива при т< I (1 - 4a)— I . При
I h21
l/4<a<l/2 схема неявная, безусловно устойчива, у=0(й2+т2), z=0(h2+т2).
Исходную схему (14) можно переписать в виде
(E-т2a(Л1 +Л —))yk+1 = [2E + т2(1 -2a)^ +Л—)]yk - (1б)
— k1 , (1б) -(E-т a(Al + Л2))y +т ф.
Оператор B=E-т2a(Л1+Л2) можно приближенно заменить факторизованным оператором С
C = (E-т 2a^)(E-т ^Л 2) = E-т 2a (Л1 +Л 2) + т 4aЛ1Л 2) = B + 0(т4) . (17)
Заменяя в (1б) оператор В на С, получим факторизованную схему
(E - т^Л1)(E - тЪЛ2)ykdl = [2E + т2(1 - 2a) (Л1 + Л2)]yk - (E - т 2a(^ + Л2))yk-l + т 2ф . (18)
Преобразуя факторизованную схему (18) к форме типа (14) и учитывая соотношение (17), получим
уь = (Лі +Л—Xay^1 + (1 - 2a) yk +ayk-1) + ф-т 2aЛlЛ -у^ ,
что отличается от схемы (14) на члены 0(т2). Поскольку схема (14)-(15) имеет второй порядок аппроксимации, то факторизованная схема (18) также имеет аппроксимацию 0(т2+й2). Таким образом, из сказанного выше следует безусловная сходимость факторизованной схемы (18) со скоростью у=0(й2+т2) при l/4<a</2.
Решение разностных уравнений сводится к последовательности одномерных прогонок по направлениям х1 и х2. В самом деле, факторизованный оператор С есть произведение одномерных трехточечных операторов E-т^Л^ а каждый такой оператор обращается одномерной прогонкой, тем самым схема (18) экономична.
Так же, как в случае одномерного волнового уравнения, задачу (12)-(13) можно решить методом прямых с применением к системе дифференциально-разностных уравнений (количество уравнений равно 2(N-l)2) метода Рунге-Кутта четвертого порядка с аппроксимацией по пространственным переменным х1, х2, 0( I h4 I ) и 0( I h2 I ). Исследования на модельных задачах показали, что оптимальным является метод прямых с аппроксимацией по пространственным координатам 0( I h4 I ).
Рассмотрим следующие модельные задачи, соответствующие точным решениям:
1а U(х1,х2,-) = (х12 + х—)-2 + 5(х1 + х2) +1, D = {0 < х1,х2 < 1; 0 < t < 2},
2а U(х1,х2,-) = ехр(х1,х2,-), D = {0 < х1,х2 < 1; 0 < t < 2}.
При к=0,1 для аппроксимации 0( | к4 | ) число уравнений 2(#-1)2=2-92=162. При к=0,05 для аппроксимации 0( | к2 | ) число уравнений 2(#-1)2=2-192=723, т.е число уравнений при аппроксимации 0( | к2 | ) в 4,5 раз больше, чем для 0( | к4 | ).
Из табл. 3 на основании проведенных исследований для волновых уравнений можно сделать вывод, что оптимальным является метод прямых с методом Рунге-Кутта 4-го порядка с аппроксимацией у=0(к4+т4).
Таблица 3
Влияние порядка аппроксимации на точность метода
Метод одном* у=0 эрных прогонок Ь2+т2) Мете \ >д Рунге-К |/=0(Ь2+т4 т т ш Метод Рунге-Кутта у=0(Ь4+т4)
Ь т И к,т И к2т Ь т 1И1(2) Ь т 1И1(2)
0.01 0.01 0.05 0.01 0.1 0.001 0.002 0.1 0.001 0.0001
0.005 0.001 0.02 0.005 0.05 0.0005 0.0001 0.1 0.0005 0.00004
3. Начально-краевая задача для уравнения колебания пластинки:
д2и д4и 0 д4и д4и
+ ^ 2~, 2 + У (X1, Х2, ..„ч
дt д х д х1 д х2 д х2 (19)
Б = {0 < х1,х2 < 1, 0 < t < Т} = Б + Г.
Начальные условия
и(х,0) = и0(х), ди(х,0) = и1(х), 0 < х = (х1,х2) < 1 . (20)
д^
Краевые условия двух видов (п - нормаль к Г):
д 2и
и Г =
и
дп2
д и
= 0 , (21)
Г
дп
=0 . (22)
В предыдущих задачах сравнивались безусловно устойчивые неявные схемы с методом прямых решения дифференциально-разностных уравнений с использованием явного метода Рунге-Кутта 4-го порядка. Теперь сравним явные разностные схемы с аппроксимацией
2242 2444
у=0(к +т ); у=0(к +т ) с методом прямых с погрешностью у=0(к +т ); у=0(к +т ).
Для аппроксимации на сетке ш'йт = {хгу = (х[г), х21)), tk = кт, /, 1 = 0, N, к = 0,М} частных
производных по пространственным переменным, входящим в уравнение (19), потребуется введение законтурных узлов, для схемы у=0(к2+т2) - один ряд законтурных узлов, для схемы у=0(к +т ) - два ряда законтурных узлов, которые определяем из разностных аппроксимаций соответствующих краевых условий (21), (22). В случае двух законтурных узлов используем два вида разностных аппроксимаций: на симметричном и несимметричном шаблонах (5), (6). В результате аппроксимации на расширенной сетке
юйг = х,*, 1 = 0, п, п = N +1; tk = кт, к = 0,М} (1=2 для схемы 0(к )
Г
(/=4 для схемы 0(к4)), получим систему явных разностных уравнений относительно сеточной
функции y..+1 = y(xf), x21), tk+1 ) .
yk+1 - 2yk + yk 1 т2
= -(yk + 2Л"1,2yk + Л«іyk)+ q , (23)
где Л^, Л^2, Лх2 - разностные операторы, аппроксимирующие частные производные в уравнении (19) с порядком 0(к2) (/=2) либо 0(к4) (/=4). Начальные условия аппроксимируем с порядком 0(т2)
у° == и0( х«,х21)) = и01
y1 - y0 т / 4 ^ 2 _ (24)
= U (i,j) + 2 +(i, j,0) - V4U0)+ O^2 ),(i, j) = 1, N -1.
Краевые условия (21), (22) аппроксимируем на расширенной сетке ш,+ с порядком
O( \ h \ ) при /=2 и с порядком O( \ h4 \ ) при /=4. Так как схема (23) и (24) явная, условно устойчивая, выбор шага т по переменной t осуществляем по правилу Рунге на последовательности сгущающихся сеток (т,т/2,т/4,...); добиваемся совпадения решений с требуемой точностью на двух последовательных сетках. Сравнение явных разностных схем с порядком O( \ h2 \ +т2) и O( \ h2 \ +т2) проводим на модельных задачах, соответствующих точным решениям:
1б U(x1,x2,t) = cosп t sinпx1 sin пx2;]
f для краевого условия (21)
2б U(x1,x2,t) = (1 +t)sinпx1 sinпx2; J
3б U(x1,x2,t) = cosпt sin2пx1 sin2nx2;]
f для краевого условия (22)
4б U(x1,x2,t) = (1 +t)2 sin2пx1 sin2пx2. J
Применяя метод прямых к задаче (19)-( 21), получим систему дифференциально- разностных уравнений, которую решаем методом Рунге-Кутта 4-го порядка, при этом используем разностные операторы по пространственным переменным порядка O( \ h2 \ ) и O( \ h4 \ ). Здесь выбор шага т также проводим на последовательности сгущающихся сеток (т,т/2,т/4,.). Сравнение метода Рунге-Кутта с погрешностью аппроксимации O( \ h2 \+т4) и
O( \ h4 \ +т4) проводим на модельных задачах 1б-4б.
Результаты исследований приведены в табл. 4, 5.
Таблица 4
Влияние порядка аппроксимации на выбор шагов h и т в явной схеме
№ задачи Явная схема
0(^+т2) 0(^+т4)
h N т M h N т M
1б 0.02 50 0.00005 20000 0.1 10 0.0012 B00
2б 0.02 50 0.00005 20000 0.1 10 0.0006 1600
Зб 0.017 60 0.000035 2BB00 0.05 20 0.00015 6400
4б 0.0125 B0 0.00003 32400 0.05 20 0.00015 6400
Таблица 5
Влияние порядка аппроксимации на выбор шагов Ь и т в методе Рунге-Кутта
№ задачи Рунге-Кутта
0(Ь2+т4) 0(Ь4+т4)
Ь N т М Ь N т М
1б 0.02 50 0.0001 10000 0.1 10 0.0025 400
2б 0.02 50 0.0001 10000 0.1 10 0.0025 400
3б 0.017 60 0.00007 14400 0.05 20 0.0003 3200
4б 0.014 70 0.00005 19600 0.05 20 0.0003 3200
Здесь погрешность на точном решении задается равной 0.001 при Т=1. На основании проведенных исследований [7] сделан вывод, что оптимальным алгоритмом является метод Рунге-Кутта с невязкой у=0(И +т ), т.к. он позволяет брать более крупную сетку по И по сравнению с 0(И2+т4), что приводит к уменьшению порядка системы в 4,5 раза. Этот алгоритм был применен для решения более сложной задачи.
4. Уравнения движения гибких пластинок в кинематической гипотезе Кирхгофа:
д2 w дм
----2—+ ^-------
дл2 дл
1
12(1 -V2)
(25)
1
V ^2 Р = - - Ь(у, м).
В уравнении (25) (х, у) е G = {0 < х < 1,0 < у < 1}, Х = 1,0 < I < 1п, Дм,Р) - известный нелинейный оператор
т, ^ д2м д2Р _ д2м д2Р д2м д2Р Ь(м, Р ) =— ------ 2-------------+ -
дх ду дхду дхду ду дх
Начальные условия:
w\t=0 =Ф:( x, у),
дм
д
= Ф 2(X, у).
(26)
Выпишем для области {0<х, у<1} один из вариантов краевых условий, при х=0; 1 (хоу). Например, шарнирное опирание на гибкие нерастяжимые в касательной плоскости ребра:
дх
2
дх
2
(27)
Здесь выбор шага И по пространственным координатам и шага т по времени в методе Рунге-Кутта осуществляется по правилу Рунге, т.е. добиваемся совпадения решения на сетке с шагом И и И/2, и соответственно с шагом т и т/2. В результате этих экспериментов был выбран шаг И=1/8 и т=0,0002 в схеме повышенной точности.
5. Слабонелинейное многомерное уравнение теплопроводности
В [1] даны схемы второго и четвертого порядков точности для стационарного и нестационарного линейных уравнений теплопроводности. В [2] приведена схема второго порядка для слабонелинейного стационарного уравнения теплопроводности в прямоугольной области в случае первой краевой задачи. Доказана сходимость к решению дифференциальной задачи. В случае третьей краевой задачи методом аппроксимации квадратичного функционала в [3] построена
г =0
разностная схема повышенной точности для стационарного линейного трехмерного уравнения теплопроводности с постоянными коэффициентами при старших и младших производных, показана сходимость к точному решению и эффективность на модельных задачах. Поэтому в данной работе мы ограничились рассмотрением первой краевой задачи. В [4] построены разностные схемы повышенной точности для многомерных линейных стационарных уравнений теплопроводности со смешанными производными, с переменными коэффициентами при младших и старших производных. В [5,6] построены разностные схемы повышенной точности для слабонелинейного стационарного многомерного уравнения теплопроводности и найден экономичный алгоритм решения систем нелинейных разностных уравнений.
В данной работе построены неявные схемы с весами и схема повышенной точности для многомерного слабонелинейного уравнения теплопроводности, на модельных задачах выявлен экономичный алгоритм решения соответствующих нелинейных разностных уравнений - это метод переменных направлений для схемы повышенной точности.
В области О = {0 < ха < 1а, а = 1,2,3} и 0<г<Т рассмотрим слабонелинейное уравнение теплопроводности
ди/дг = £ д 2и/ дх2а - Ко {х,г,и,ди/дх), х = (, Х2,Хз )е О, г > 0,
а=1 (28)
ди/дх = (ди/дх1, ди/дх2, ди/дх3),
с краевыми условиями первого рода и(х,г)| = д (х,г), Г - граница области О, и начальным
условием и(х,0)=и0(х) х е О .
Слабонелинейность уравнения (28) означает [2], что функция К0(х, г, р0, р1, р2, р3), где р0=и, ра=ди/дха, а=1,2,3, определена при х е О, г е [0,Т], |р0 I, |р11, |р2 I, |р31 и непрерывна по х, г при фиксированных ра(а=0,1,2,3), а также существуют производные от функции К0 по ра, которые удовлетворяют условиям
М0 > дК0/др0 > 0, |дК0/дра\ <М, а = 1,2,3 .
5.1. Явная разностная схема. Введем прямоугольную равномерную сетку
®И = { *2 К 13И3 )} } = 0, Ма , Иа = 1а/ Ма , а = 1,2,3 и ®т = {г„ = nт}, П = 0,т,Т = Т1™ и
сеточную функцию у = уп = у ((, /2И2, /3И3, гп). Перейдем к безындексным обозначениям, полагая У = у”+1, уг = (у?- у)/т . Нижний индекс будем указывать только в том случае, если он отличен от /а (а=1,2,3), например, уха = (у -уг-а-1)/Иа. Рассмотрим разностные операторы
23
ЛаУ = Ухаха = (Уга-1 - 2У + У,а+1)/К , Лу = £ЛаУ . Поставим в соответствие задаче (28) разно-
а =1
стную схему
У{ = ЛУ - 0.5[K0 (Х, г, У, Ух ) + K0 (Х, г, У, Ух )] Ух = (Ух1 ,Ух2 ,Ух3) . (29)
Краевые и начальные условия выполняются точно. Непосредственно можно убедиться, что погрешность аппроксимации разностной схемы
у = Ли-КДх,г,и,их)-и{ = О (\И2\+т) ,
где
К (х, г, и, их) = 0,5 [к 0 (х, г, и, их)+K0 (х, г, и, их)]. (30)
Разностная схема (29) является явной и условно устойчивой (т<И2/6, И=шт Иа), требует малых шагов И и т, поэтому применяется редко [1].
5.2. Схемы с весами. Рассмотрим неявную схему с весами для уравнения (28)
Уг =л( +(1 -о)у)+ф, х еШи, г е®т, (31)
ф = -К1 = -К1(хг”y,ухX г = гп + V2, у = у(x,0.
Для оценки порядка аппроксимации схемы (31) представим о? + (1 -о) и = (и+и )/2 + (о-0,5) т и,, иг = и + О(т2), Л„£7 = Ьр + И’/12Ц и + О (И4)
и вычислим невязку
у = Л(о и? + (1 -о) и) + ф- и{ =Л (г? + и )/2 + (о-0,5) тЛиг +ф-иг =
= (и -К1 -й7)+ (о-0,5)тШ + (ф + К)+О(|И|2 +т2) =
= (о-0,5)тЬи +(ф + К1)+О(|И|2 +т2), где Ьи = Ли.
Таким образом, у=О( | И | 2+т2) при о=0,5, ф = -К1; у=О( | И | 2+т) при о^0,5, ф = -К1. Вместо схемы (31) можно рассматривать схемы с различными весами оа по направлениям ха
У( =1 Ла (оаУ + (1 -оа )У)+ф, ф = -К1. (32)
а=1
5.3. Схема повышенной точности. Построим схему повышенной точности для уравнения (28). Покажем, что схема
3 3 1^3
У( =ЕЛа(°а.£ + (1 -оа )у)+Е Иа/12 ЕЛаЛрУ +ф , (33)
а=1 а=1 р^а
где
оа= 0.5 - И;/(12т), ф = -ГК +1 И;/12Л„К,) , (34)
имеет погрешность аппроксимации у=О( | И | 4+т2).
Запишем невязку
¥ = £( (( + и)/2 + (о„- 0,5)тЛи и,)+ £И;/) ЛаЛри и,.
а=1 а=1 р^а
Учитывая, что
($+и)/2 = и+о(т2), лаиг = ьаи+о(т2 + и2), Ци = ци -ь1ь2и -ьци + ь1К1, ь\и = ь2и - ь2 ь1и - ь2 ь3и + ь2 К1, ь2ъи = ь3и/ - ь3 ь1и - ь3 ь2П + ь3 К1,
и/ - ьи + К = 0,
получим
у = £ [и[12 + К- 0,5)т]ьа Г7 + ('ф + К + £ И;/12 Ц К') + О (т2 + |И|4).
а=1 V а=1 )
Таким образом, у=О(т2+ | И | 4), если оа и ф взять согласно (34), при этом функция К0 в уравнении (28) должна иметь непрерывные производные по ха до четвертого порядка. Для вычисления Ла К1 в (7), (30) используем формулы
Ла к=(к о; -1)- 2К+к (/„+1)^ иа,
где
*,(4-1) = К ((і а - 1*, іп, и ((і.- 1)Ла), и, (/„-1)).
Исследования на модельных задачах показали, что такая разностная производная от нелинейной функции имеет второй порядок точности, что и требуется для схемы (33).
5.4. Экономичные алгоритмы. Для решения систем разностных уравнений (32), (33) на каждом временном слое іп можно использовать метод матричной прогонки. Однако, он требует большого числа действий 0(п4()), где п0 - число узлов сетки ш* (по=#з, где
#=Лга, а=1,2,3). Поэтому непосредственное использование схем нецелесообразно. Каждая из схем (32), (33) может быть заменена схемой того же порядка аппроксимации, но требующей для определения У последовательного применения скалярной прогонки для трехточечного уравнения и затраты О(п0) арифметических действий. Такая схема называется экономичной - это метод переменных направлений. При ее написании используется схема повышенной точности (33). Приведем схему переменных направлений для двумерного уравнения теплопроводности [2]:
2—^ = аДУ +(1 -о2)у + , = (1 -а!)Лу + а2Л2У +(1 -а1)фп,
т т
где У = УЫпX У = У(хЛ +т/2), у = У(х , гп Егат.
Запишем краевые и начальные условия
У(х,0) = ио( х), X Ей*, У = Д ^) при і2 = 0М2, У = Ц при і1 = 0,Nl,
где
Ц = а1 Цп+1 + (1 -а1) Дп -тЛ2(а1 а2 Цп+1 - (1 -а1)(1 -а2)Цп ), а. = 0,5-*=/(12т), ф" = -(*1 + *2/12ЛК + *2/12Л,^).
Приведем одну из модельных задач, на которых проводились исследования.
и (х1,х2,ї) = Ї 2[( х1 + 0,5)2 + (х2 + 0,5)2] +1 - точное решение задачи (28) в двумерной области О = {0 < ха < 1, а = 1,2}, ґє[0,1]; этому решению соответствует
К0 = -2 [(х1 + 0,5)2 + (х2 + 0,5)2] +16(р - 1) ^Др2 + Р2 ) , где р=и, ра=ди/дха. Для достижения точности 10-4 потребовалось для явной схемы в 370 раз больше действий, чем для схемы повышенного порядка точности, а для неявной схемы (а1=а2=1) в 9 раз больше, чем для схемы повышенного порядка точности.
Вычислительные эксперименты на модельных задачах, для которых известны точные решения, показали, что самым экономичным алгоритмом является метод переменных направлений для схемы повышенной точности, т. к. позволяет использовать более крупную сетку и сокращает порядок системы на каждом временном слое в десятки раз по сравнению со схемой второго порядка, что очень важно при решении многомерных задач. Разработанные алгоритмы и комплекс программ представляет интерес для практики решения нестационарного слабонелинейного многомерного уравнения теплопроводности в прямоугольной области, т.к. требуют только задания функции К0, границ прямоугольной области и числа разбиений по пространственным и временной координатам.
Итак, на основании исследования целого ряда задач для математической физики можно сделать вывод, что использование разностных схем повышенной точности по пространственным координатам приводит к значительному сокращению порядка систем разностных уравнений на каждом временном слое, что особенно важно при решении многомерных задач.
ЛИТЕРАТУРА
1. Самарский А.А. Теория разностных схем. М.: Наука, 1977. 580 с.
2. Самарский А. А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 244 с.
3. Вахлаева Л.Ф., Крысько В. А .Устойчивость гибких пологих оболочек в температурном поле // Прикладная механика. 1983. Т. 19. № 1. С. 16-23.
4. Вахлаева Л.Ф. Разностные схемы повышенной точности для эллиптических уравнений // Математика и ее приложения: Межвуз. сб. науч. тр. Саратов: Изд-во Сарат. ун-та, 1991. Вып. 2. С. 74-77.
5. Вахлаева Л. Ф., Вахлаева Т. В. Разностные схемы повышенной точности для слабонелинейного эллиптического уравнения и сравнение со схемами второго порядка точности / Сарат. гос. техн. ун-т. Саратов, 1996. 9 с. Деп. в ВИНИТИ 24.09.96 № 2855-В96.
6. Вахлаева Л.Ф., Вахлаева Т.В., Павлова Е.А. Экономичные алгоритмы решения разностных краевых задач для эллиптических уравнений // Математика. Механика: Сб. науч. тр. Саратов: Изд-во Сарат. ун-та, 2001. Вып. 3. С. 21-23.
7. Вахлаева Л.Ф., Вахлаева Т.В., Назарьянц В.О. Исследование разностных схем второго и четвертого порядков точности для нестационарных уравнений математической физики // Математика. Механика: Сб. науч. тр. Саратов: Изд-во Сарат. ун-та, 2000. Вып. 2. С. 159-161.
Крысько Вадим Анатольевич -
доктор технических наук, Соросовский профессор,
Заслуженный деятель науки и техники РСФСР, заведующий кафедрой «Высшая математика»
Саратовского государственного технического университета
Вахлаева Людмила Федоровна -
кандидат технических наук,
доцент кафедры «Математическая физика и вычислительная математика»
Саратовского государственного университета им. Н.Г. Чернышевского
Молоденкова Татьяна Викторовна -
кандидат физико-математических наук, доцент кафедры «Высшая математика»
Саратовского государственного технического университета
