CC BY
42
УДК 519.633+519.635+631.41
Вестник СПбГУ. Сер. 1. Т. 3(61). 2016. Вып. 2
ОБ ОДНОМ МЕТОДЕ ЧИСЛЕННОГО РЕШЕНИЯ СИСТЕМЫ ПАРАБОЛИЧЕСКИХ УРАВНЕНИЙ*
Н. Н. Семенова1'2, В. В. Терлеев3, Г. И. Сухорученко2, Е. Е. Орлова1, Н. Е. Орлова1
1 Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9
2 Всероссийский исследовательский институт защиты растений (ВИЗР), Российская Федерация, 196608, Пушкин, шоссе Подбельского, 3
3 Санкт-Петербургский политехнический университет Петра Великого, Российская Федерация, 195251, Санкт-Петербург, Политехническая ул., 29
В работе приводится алгоритм численного решения системы двух параболических уравнений специального вида и доказывается его устойчивость. Даются оценки полученных численных решений. Определяются условия, гарантирующие применимость процедур прогонки и дается оценка на шаг дискретизации по времени, при котором итерационный процесс имеет внутреннюю сходимость. Приводится пример использования алгоритма в почвоведении. Биб-лиогр. 10 назв. Ил. 2.
Ключевые слова: метод прогонки, итерационные процедуры, весовой коэффициент, монотонность, почвоведение.
1. Введение. Рассматривается система дифференциальных уравнений в частных производных параболического типа в прямоугольной области. Коэффициенты, входящие в эту систему уравнений, зависят от искомых переменных, поэтому данная система дифференциальных уравнений является квазилинейной. Нахождение аналитических решений таких систем даже в прямоугольной области возможно лишь в некоторых частных случаях, поэтому в работе использовался традиционный численный метод конечных разностей (см., например, [1]), дополняемый итерационным процессом. Численному решению параболических уравнений уделялось достаточно много внимания [2, 3], однако практически не рассматривались системы параболических уравнений, которые получили распространение в связи с задачами, возникающими в почвоведении и геологии [4, 5]. Решение этих задач требует также высокой точности численных решений. В работе предлагается способ аппроксимации производной первого порядка, что позволяет построить схему второго порядка точности по переменным, входящим в данную систему, и получить оценки решений системы, а также разрабатывается метод матричной прогонки. Расчетный профиль разбивался на N слоев с толщиной Н¿, причем Н не обязательно одинаковы, а временной промежуток— на периоды длительностью т. Тогда систему дифференциальных уравнений возможно аппроксимировать дискретным аналогом по конечно-разностной схеме с весами [3].
2. Основные понятия и обозначения. Построение расчетной схемы. Дана система двух параболических уравнений специального вида с неизвестными функциями Р(х, £), Ст(х,Ь) от пространственной переменной х и временной переменной £
(1)
* Работа выполнена при финансовой поддержке РФФИ (грант № 16-04-01473-а). (¡5 Санкт-Петербургский государственный университет, 2016
Здесь
яд (яр \
ст = сь(в + Рк-1), КР) = ор, ^ = к(р,сь)(—-1у
/п = — KdegCт = —Хь Сь — Хб СБ ,
где /б является суммой функций: функции /р, представляющей процесс кинетики первого порядка переменной Ст(х,Ь) с относительной скоростью К^, и непрерывно дифференцируемых функций /е и /р (сут-1). Считается также, что функции (0), ц(Р), К(Р,Ст), ) удовлетворяют условию Липшица по переменным Р и
Ст е некоторой константой Ь (е учетом существования 0(Р) и 0-1(Р)). Тогда можно считать, что все функциональные зависимости, входящие в уравнения системы (1) удовлетворяют условию Липшица с константой Ь > Ь.
Введем следующие обозначения: I — длина расчетного слоя по пространственной переменной х; (Ь) и (Ь) —потоки по переменной Р, см-сут-1, а ¿ор(Ь) и ¿ьр (Ь) —потоки по переменной С соответственно через верхнюю и нижнюю границы, мг-см-2-сут-1; С01,г —дискретное распределение начального содержания по переменной Сь, мг-см-3; —дискретное распределение начального содержания для потока Лр, см-сут-1, ¥Ггг —доля потока, для верхней границы расчетного слоя (Рг,г = 1 при некоторых условиях); Ьг —моменты задания дискретного распределения, сут; г — номер задания (г = 1,...,п, где п — число заданий); 5(Ь) —импульсная функция (¿-функция Дирака). В данном случае использование функции 6(Ь) позволяет записать эффект точечного (дискретного, в моменты времени Ьг) задания переменной Сь (считается, что в моменты времени Ьг Ст = Сь, т.е. Сб(Ьг) = 0). Тогда верхние граничные условия имеют вид
7^1 _ т,(дР
■7 |*=0 = ~К ~ 1) = ^М*) + Ъ ~
Мх=о = ^нь ^ + = -ТОР®' (2)
¿ор (Ь) = ^^ Рг,г С'ь,г$(Ь — Ьг),
г=1
а нижние граничные условия определяются как >х-1 \дх
(3)
М^ = ВН, (в, -п + |а=, = ыь).
Будем считать Ьо(г-о) начальным моментом расчета, а функции (Ь), (Ь), ¿ьр (Ь) известными. Предположим также, что начальное распределение переменной Р по пространственной переменной х задается функцией Ро(х). Будем также считать, что в начальный момент времени Ьо функция Со (х) (0 < х < I), описывающая распределение переменной Ст по пространственной переменной х, может быть отлична от нуля только в некотором ограниченном интервале, т. е. будет задаваться ступенчатой
функцией вида
\Сы, 0 < х < Ьо,г = 0; Ьо > 0, Со (х) = <
\0, х > 0, (4)
Со (х) = Сы6(х), Ьо = 0, I = 0,
где Сыо, Ьо —константы и возможен случай обращения в ноль величины Ьо, а 6(х) — по-прежнему импульсная функция (¿-функция Дирака).
Для численного решения системы (1) с граничными условиями (2), (3) и начальными условиями (4) воспользуемся конечно-разностными методами и предложенной модификацией консервативных разностных схем с переменным шагом по пространственной координате [6].
Введем в рассмотрение следующие сетки по пространственной переменной х: с^о = {х1,..., ждг} с шагами К\ = Щ-, Ъ^ = — х^, /гдг = хт~хм ; = {хо = х\/2 = 0; ^¿+1/2 = х1_1/2 + К] хт = жлг+1/2 = 1]Ы = ипо времени и шТ = {Ь = ]т, т >
0, ] = 0,1,...,К; К = [Т/т]}, где I — длина расчетного слоя, Т — рассматриваемый промежуток времени [1, 2]. Граничные условия третьего рода (2), (3) выбраны как наиболее общий тип граничных условий.
В целочисленных узлах (х^,1?,1 = 1,...,^ ] = 1,...,К) вычислим сеточные функции Р?, С]ы, соответствующие функциям Р(х,Ь) и Сы (х,Ь), а в узлах (х^+1/2) —сеточные функции, являющиеся аналогом потоков в первом и во втором уравнениях системы (1):
X = Xi+1/2,t = tj
дС
^Jr\n = DHL{()1Jw)—±-JwCL =Jp{xi+i/2,tj), j=TJ{.
' 1 dx X = Xi+i/2,t = tj
(5)
Рассмотрим конечно-разностную аппроксимацию системы (1) схемой с весами (6 — весовой параметр, 0 < 6 < 1), которая при 6 = 0 представляет явную схему, при 6 = 1 —неявную, а при 6 = 1/2 — схему Кранка—Николсена. Соответствующая система разностных уравнений, аппроксимирующая систему (1), в матричной форме имеет вид
Р з+1 _ р з
——+ 6А{+6рз+1 = - (1 - 6)А{+6рз, + ЗА^С^1 = - (1 - 5)А?+*С3Т,
где Рз = (Р?), СТ = (С?^, I = 1,...,^ ] = 1,...,К, — ^мерные векторы, а матрицы-функции имеют вид
Ч
где
Al+S =
(Л™
A12) | , (7)
\Л(3)
= -Ж ...
\ hihi ' hihi '
где
А(2)=
А
(2)
0, -
к1
-1/2
К1
k¿—lk¿
■1/2
+К
+ 1/2
к
+ 1/2
k¿—1 k¿
k¿—lk¿ '
А
(3)
0, ...,--
А22+^
к
N-1/2
к
N-1/2 hN-lhN ' ft.JV-lft.JV
а22) 23)
А
а21) =
3/2 /11/11
+
3/2
к!
- К3+6
г,4ед'
В:
1+6 3/2
/ц/ц '
...,0
-1/2
к^-1 к
-1/2
В
1+6 1/2
2к»
к®_1 к»
+ -
0:
+ 1/2
к» к»
+ -
¿-1/2 ¿+1/2 2к»
. П1+» 7
_Т(3 + 6 ¿+_1/2 | »+_1/2 Мед> 2ft¿
; = О ^У-1/2 ^-1/2 _ ^' + <5
1 ' ftJV —lftJV, ft.jv-ift.jv ftJV г^ед
М 3+6
Ф
ь
(1=3,3 + 1)
(3+6
0 0
/я{+6
0 0
0
3+6
К
0 0
■3+6 N
(9)
(10)
Входящие в определение этих матриц (7)—(10) величины находятся из выражений
,3+6
п3'+6 = п„г (03+6 ¿УУ,3+6) + N3+° . N3 + ° г+1/2 = ^нь(дг ,-^г+1/2 ) + ^ г+1/2. г+1/2
30
— I — О:>+Й 4- К ап - К] + ° — К( Р-
Яр \р=р'+6 . Кг = г + Кг,Брг. Кг+1/2 = К (Рг
3+6
3+6 ) +1/2)
Л3 + 6 т™,3 + 6
3+6
3+6
Л.
ж,3+6 Л +6
г+1/2
803+6 03+6
+1/2
Р3+6 = Р3+15 + Р3 (1 — 6). 03+6 = 0(Р3+6). СТ+г6 = СТ+/5 + 3г(1 — 6);
4+1/2 21г+1 ' 2К+1
Кг,+6д = Хь 03+6 + Хб Кг Брг.
Вектор-функции Р3+6, 03+6 записываются как
к
р з+6
к3
3 + 6 _ т Ш,з + 5
3/2
Тч
■1/2
К3
■1/2
V
к1 + 6 + т Ш,1 + 6 N—1/21
клг
С3+6
/
(11)
(12)
0
1
1
о
7
0
к
к
При записи системы (6) использовалась конечно-разностная аппроксимация второго порядка точности по х граничных условий, задаваемых уравнениями (2), (3) на четырехточечных шаблонах, состоящих из узлов {(х1,Ь?+1), (х2,Ь?+1), (х1,Ь?), (х2,Ь?)} для верхней границы и из узлов {(ху-1,Ь?+1), (ху,Ь?+1), (ху-1,Ь?), (ху)} для нижней границы. Будем исследовать разностную схему (6) с соответствующими граничными условиями и весовым параметром 6 как наиболее общую задачу.
Система (6) является системой квазилинейных разностных уравнений, для решения которых прибегают к итерационному процессу: на каждом шаге разностная схема оказывается линейной и становятся применимыми методы линейной алгебры. В качестве начальной итерации берется решение, найденное на предыдущем временном слое. Применение схемы Кранка—Николсена (при 6 = 1/2) в сочетании с соответствующей аппроксимацией первой производной дбС/ дЬ дает схему второго порядка точности по переменной Ь [7]. При 6 =1 получаем неявную схему первого порядка точности по Ь.
Очевидно, первое разностное уравнение системы (6) аппроксимирует первое уравнение системы (1) со вторым порядком точности по х. Для того чтобы получить аппроксимацию второго порядка точности по х и для второго разностного уравнения системы (1), производная Сы/дх была аппроксимирована по формуле
д^Сь !)
дх ~ 2Ы 2йг~1
Однако такая схема не является монотонной, и для того чтобы применить метод прогонки для решения системы (1), требуется дополнительное исследование. Рассмотрим еще вариант конечно-разностной аппроксимации Сы/дх [1]
дх ^ +
Аппроксимация выражения (14) содержит односторонние производные, учитывающие знак функции Л^, и дает устойчивую монотонную схему, для которой всегда применимы формулы прогонки, однако такая схема имеет первый порядок точности по х. Расчет по такой схеме может привести к значительным искажениям численного решения системы (1).
Линеаризация системы разностных уравнений (6) для осуществления итерационного процесса приводит к системе вида
гг3+1,з + 1 _ из
В{рз+*'а)-+ <ШР+1'8+1 = г{Р*+*'а) - (1 - 5)Аи*,
т
и =(фСт) , * =( с), ^ = и3, (15)
(И (Р) 0 N Л_(А1 0
В(Р) = X ' Р - А
V 0 Е1' 10 А
2
Оставив пока вопрос о сходимости итерационного процесса открытым, исследуем устойчивость соответствующей (15) разностной схемы. Данная схема является линейной относительно неизвестных векторов Р3+1,3+1, С3+1'>1+1, и, не умаляя общности, можно изучать ее устойчивость при в = 0.
Решим систему (15) относительно неизвестных векторов р3+1,1 = Р3+\ СТ+1,1 = С 3+1:
Ст
Р3+1 = Т2Р3 + т (Е + 6тА31М3 )-1Р3, С3+1 = Т1С3 + т (Е + 6т (Ф3+1 )-1А2 Ф3+1)-1С3,
где Т1,Т2 — матрицы перехода
(16)
Ti = (E + ST(Ф3+ )-1A2) (Ф3+ )-1 Ф3 (E - (1 - S)T(ф> )-1A2 T2 = (E + SrA{M3) 1 (e - (1 - S)rA{M3) ,
(17)
Согласно работе [8] для устойчивости разностной схемы достаточными являются следующие условия:
1|Т1|| < 1, ЦТ2Ц < 1, (18)
где ||*|| — норма матрицы (здесь и далее знак нормы обозначает норму Фробениуса). Если условие (18) заменить более слабым
||Т1|| < 1 + 0(т), ||Т2|| < 1 + 0(т), (19)
определение устойчивости при малых т будет допускать экспоненциальное возрастание по времени погрешности округления. Найдем условия, при которых схема (15) будет устойчивой в смысле неравенств (18).
3. Теоремы об устойчивости разностных схем
Теорема 1. Если функция (х,Ь) монотонно убывает по х при всех Ь € [0, Т], а 0(х,Ь) —убывающая функция времени, разностная схема (15) устойчива и справедливы следующие оценки:
р 3+1
г3+1
гт
< ||PU|| + тт1,
ч0 | 'Т
< ||гТ|| + ТШ2, j = 1,...,K,
m1 = max [sup \J0W (t)\, sup \JW (t)\j , (20)
m2 = max (sup J (t)Co(t)\, sup \JW (t)C (t)\
V t t
где знак нормы обозначает норму Фробениуса.
Предположим теперь, что функция JW (x,t) произвольна. Пусть имеем N * = 2mrAXi,3 \ JSW-1Lj2\, тогда матрица A2 будет полуограниченной при всех j, и будет выполняться следующее неравенство:
A2 > -N*E. (21)
Нетрудно показать, что в этом случае вместо первого условия (18) будет выполняться первое условие (19).
Теорема 2. При выполнении неравенства (19) и в предположении убывания функции О(х,Ь) по времени разностная схема (15) устойчива в смысле неравенств (19) и выполняются следующие оценки:
\Р3+1 II < ||Р°|| + тт1,
С3+1 Ст
, 2К* < ехр ( -=-
( М А
в = тах —---—
з \ 1 + 0.5г /
С°
ехр
то2
1,...,К,
(22)
где т1 и т2 те же, что и в условиях (20).
Для того чтобы убрать в условиях теорем 1 и 2 ограничение на функцию О(х,Ь), вместо разностной схемы (15) рассмотрим схему, в которой матрица А и вектор * останутся прежними, а в качестве матрицы В и вектора V возьмем следующие матрицу и вектор:
Р \ п _ (И (Р3+6) 0 \ С?); В V 0 &+ч;
и
(23)
*
О
Очевидно, что для разностной схемы (23) характер убывания или возрастания О(х, Ь) на устойчивость не влияет.
После того, как установлены условия устойчивости схемы (15) и произведены оценки вычислительной погрешности при переходе с одного временного слоя на другой, покажем корректность вычислительного процесса на данном временном слое. Не умаляя общности, можно исследовать этот вопрос также для в = 0.
Будем решать систему (15) методом прогонки последовательно. Применимость формул прогонки для определения вектора Р3+1 достаточно очевидна. Из положительности матрицы А31 следуют монотонность разностной схемы и применимость метода прогонки для любых Н и т.
Из анализа второго уравнения системы следует, что коэффициенты при неизвестных СТ 1 ,С? 1,С? 1+1 могут менять знак (соответствующая разностная схема немонотонна). Модифицируя теорему, приведенную в работе [8, с. 64, 65], получим результат, гарантирующий условие применимости метода прогонки и в этом случае.
Теорема 3. Если на каждом временном слое ] выполняется неравенство:
г < тах к^v з = 1,.. „К, (24)
1 \ "1 — 1/21 )
система (16) при 6 > 0 имеет единственное решение С? ^, г = 1,...,М, причем для нахождения этого решения применимы формулы прогонки.
Отметим, что схему с весами целесообразно использовать только при 6 > 1/2. Исследуем сходимость итерационного процесса, заданного матричным уравнением (16). Покажем, что при выполнении условий теоремы 1 имеется внутренняя сходимость итерационного процесса. Записав систему (16) для в-й и (в + 1)-й итераций и произведя вычитание одной системы из другой, получим последовательно следующие
неравенства:
3,а Гр
3,а
Гс
< тЬ
Е + (А3,-1)-1
<т1\\вС3'8-1\
' а-1
3,8-1
Е + (АГ-1)-1
3,8-1
гра = Р3а — Р3
Ь=—, \\М3\\>1л3, ¡3
г3са = (вСт)3,8 — (ОС)3,8-1,
где Ь — константа Липшица, определенная выше.
Из неравенств (25) следует, что всегда можно выбрать такой шаг, при котором итерационный процесс будет иметь внутреннюю сходимость, причем при выполнении условий устойчивости в смысле неравенств (19) шаг т надо брать более мелким, чем при выполнении условий устойчивости в смысле неравенств (18).
Отметим, что при значениях функции ¡(Р) близких к 0, предлагаемый итерационный процесс может давать очень большие погрешности (в неравенстве (24) величина ¡л3 может быть очень малой). Поэтому для нахождения решений первого уравнения системы (1) в этих случаях использовался метод потоковой прогонки [9].
4. Пример использования алгоритма. Чтобы оценить численную схему, заданную уравнениями (1)—(5), и соответствующую программу для ЭВМ, сравним результат численного эксперимента с данными натурных опытов по изучению перераспределения влаги и ионов хлора в метровой почвенной толще в процессе испарения с поверхности почвы [10]. Численное решение находим по конечно-разностной схеме с использованием в первом варианте аппроксимационного уравнения (13), а во втором— выражения (14). Граничные условия получим, проинтерполировав известные из полевого эксперимента значения влажности почвы и концентрации ионов хлора у поверхности исследуемой толщи и на глубине 1 м в моменты времени 3, 7, 10, 40, 100 сут. Всего рассмотрим 15 слоев по вертикали со следующими толщинами: Н = 2 см при г = 1,...,5; Н^ = 5 см при г = 6, 7; Н = 10 см при г = 8,..., 15.
Влагопроводность и дифференциальную влагоемкость почвы (т. е. функции в(Р), ¡(Р), РНЬ(в^№) и К(Р,Ст)) параметризуем следующими уравнениями:
К (Р,Ст ) = Кт(Ст) Кт(Ст) =
Р Р т
~Т
Р
п(Ст) = по + а1Ст/2 + а2Ст + аз СГ,
Риь(в, J^) = Ро(Ст)ваи + КЛед(в, Ст)^W |, 0 < КЛес(в, Ст) < 1,
где п(Ст) —динамическая вязкость раствора, Кт(Ст), Кт —коэффициенты фильтрации почвенного раствора и воды соответственно, а по —динамическая вязкость воды; причем для данного типа почвы будем считать, что по ~ п(Ст), т.е. Кт ~ Кт(Ст) = 3 см/сут и зададим следующие параметры: вт = 0.005; Рт = —2 см;
По
, в
0.9; с! = 10; Ро = 0.0002 (типичные значения для иона хлора),
Х(в,Ст) = 0.55. Тогда все условия для устойчивости разностной схемы выполнены.
с
е
в
1
а
Для внутренней сходимости метода итераций и применимости формул прогонки оказалось возможным выбрать достаточно грубый шаг по времени — т = 0.5 еут.
слой почвы, см
Рис. 1. Содержание хлора (мг/см3) в слоях почвы на третьи сутки эксперимента.
200 — "Я 180 ^
.О :
2 160 г
«f
& :
а но г
о
К .
S 120 Г 100 Г 80 г 60
& §
О
-о- эксперимент о расчет по вар. 1 -с- расчет по вар. 2
i i i i i i i i
i i i i
0-2
4-6
8-10 15-20 30-40 50-60 70-80 90-100
слой почвы, см
Рис. 2. Содержание хлора (мг/см3) в слоях почвы на сотые сутки эксперимента.
Сравнение расчетных и экспериментальных кривых на третьи (рис. 1) и сотые сутки (рис. 2) после начала опыта показывает, что расчеты, произведенные по первому варианту конечно-разностной процедуры, дают хорошее совпадение, в то время как расчеты по второму варианту имеют значительные расхождения с экспериментальным материалом.
5. Выводы. Таким образом, использованные в работе расчетные схемы для нахождения численных решений системы (1) при сделанных предположениях относительно функций, входящих в вышеуказанную систему, устойчивы. Определены условия, гарантирующие применимость процедур прогонки, и дана оценка на шаг дискретизации по времени, при котором итерационный процесс имеет внутреннюю сходимость.
Литература
1. Самарский А. А., Гулин А. В. Численные методы. М.: Наука, 1989. 432 с.
2. Бахвалов Н. С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Лаборатория Базовых Знаний, 2002. 636 с.
3. Годунов С. К., Рябенький В. С. Разностные схемы. Введение в теорию. М.: Наука, 1977. 440 с.
4. Bresler E., McNeal B.L., Carter D.L. Saline and sodic soils. Principles Dynamics Modeling. Berlin; Heidelberg; New-York: Springer-Verlag, 1982. 212 p.
5. Рогачев М.В., Примак П.Н. О решении обратной электродинамической задачи видеоимпульсного георадара // Разведка и охрана недр. 2004. № 12. С. 10—16.
6. Заславский Б. Г., Семенова Н.Н., Терлеев В. В. Численный расчет миграции влаги и растворенного вещества в почве // Физика и физико-химия корнеобитаемого слоя почвы: Сб. науч. тр. / отв. ред. О. Г. Усьяров. Л.: АФИ, 1989. С. 33-48.
7. Березин И. С., Жидков Н. П. Методы вычислений. М.: Наука, 1966. Т. 1. 464 с.
8. Марчук Г. И. Методы вычислительной математики. М.: Наука, 1977.
9. Дегтярев Л. М., Фаворский А. П. Потоковый вариант метода прогонки для разностных задач с сильно меняющимися коэффициентами // Журн. вычисл. матем. и матем. физ. 1969. Т. 9, №1. С. 211—218.
10. Семихненко П. Г., Кондратьев В. В. Особенности движения влаги в выщелоченном черноземе Кубани // Теоретические вопросы обработки почв. Л.: Гидрометеоиздат, 1973. С. 162-168.
Статья поступила в редакцию 17 сентября 2015 г. Сведения об авторах
Семенова Наталья Николаевна — кандидат физико-математических наук, доктор биологических наук, доцент; nnsemenova@yandex.ru
Терлеев Виталий Викторович —доктор сельскохозяйственных наук, профессор; vitaly_terleev@mail.ru
Сухорученко Галина Ивановна — доктор сельскохозяйственных наук, профессор; info@vizr.spb.ru
Орлова Елена Евгеньевна — кандидат сельскохозяйственных наук, старший преподаватель; orlova55@mail.ru
Орлова Наталия Евгеньевна — кандидат биологических наук, доцент; orlova55@mail.ru
ABOUT ONE ALGORITHM FOR THE NUMERICAL SOLUTION OF THE SYSTEM OF THE TWO PARABOLIC EQUATIONS OF THE SPECIAL TYPE
Natal'ja N. Semenova1,2, Vitaly V. Terleev3, Galina I. Suhoruchenko2, Elena E. Orlova1, Natalia E. Orlova1
1 St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation; nnsemenova@yandex.ru, orlova55@mail.ru
2 Institute of Plant Protection (FSBSI VIZR), Podbelskogo shosse, 3, St. Petersburg, Pushkin, 196608, Russian Federation; nnsemenova@yandex.ru, info@vizr.spb.ru
3 Peter the Great St. Petersburg Polytechnic University, Polytekhnicheskaya ul., 29, St. Petersburg, 195251, Russian Federation; vitaly_terleev@mail.ru
The work provides an algorithm for the numerical solution of a system of the two parabolic equations of the special type, and proved their sustainability. There are generated an assessment of the numerical solutions. There are defined the conditions for guaranteeing the applicability of the double-sweep method, and the estimation of the sampling step in time in which the iterative process has an internal convergence. There are made an example of how to use the algorithm in soil science. Refs 10. Figs 2.
Keywords: double-sweep method, iterative techniques, weighting coefficient, monotonicity, soil science.
References
1. Samarsky A. A., Gulin A. V., Numerical methods (Nauka, Moscow, 1989, 432 p.) [in Russian].
2. Bahvalov N.S., Hzidkov N.P., Kobel'kov G.M., Numerical methods (Laboratory Bazovih Znanij, Moscow, 2002, 636p.) [in Russian].
3. Godunov S.K., Rjaben'kij V. S., Finite-difference schemes. Introduction to theory (Nauka, Moscow, 1977, 440p.) [in Russian].
4. Bresler E., McNeal B.L., Carter D.L., Saline and sodic soils. Principles Dynamics Modeling (Springer-Verlag, Berlin-Heidelberg, New York, 1982, 212p.).
5. Rogachev M.V., Primak P. N., "About solution of the inverse electrodynamical problem of the videoimpulse georadar", Razvedka i ohrana nedr (12), 10—16 (2004) [in Russian].
6. Zaslavsky B.G., Semenova N.N., Terleev V. V., "Numerical calculation of the water and soluble substance migration in the soil", Phizika i phiziko-himija korneobitaemogo sloja pochvi, 33—48 (Ed. O.G. Us'jarov, AFI, Leningrad, 1989) [in Russian].
7. Berezin I.S., Hzidkov N.P., Computational methods 1 (Nauka, Moscow, 1966, 464p.) [in Russian].
8. Marchuk G.I., Computational mathematics (Nauka, Moscow, 1977, 456p.) [in Russian].
9. Degtyarev L. M., Favorskii A. P., "Flow variant of the sweep method for difference problems with strongly varying coefficients", USSR Computational Mathematics and Mathematical Physics 9(1), 285—294 (1969).
10. Semihnenko P. G., Kondrat'ev V. V., "Peculiarity of the water migration in the leached black (chernozem) soil", Teoreticheskie voprosi obrabotki pochv, 162—168 (Gidrometeoizdat, Leningrad, 1973) [in Russian].
