CC BY
16
УДК 519.5
ВЫЧИСЛИТЕЛЬНАЯ МОДЕЛЬ НЕСТАЦИОНАРНОГО УРАВНЕНИЯ ПЕРЕНОСА ПРИМЕСИ НА ОСНОВЕ МЕТОДА ВЗВЕШЕННОЙ
НЕВЯЗКИ
О 2004 г. В.И. Наац
The construction of calculation model for the unsteadily-state equation of source transferring in the frontier layer of the atmosphere, which is based on weigh discrepancy method, method of finite elements and approximation, minimization technique of discrepancy function, to carry out. The analytical model of calculation algorithm are receiving, some calculation property of this algorithm are discussion.
В статье излагается методика построения вычислительной модели нестационарного уравнения массопереноса на основе метода взвешенной невязки, метода конечных элементов и методов минимизации. Излагаемая методика соответствует вариационному подходу к построению математических моделей, используемых далее в задачах мониторинга и прогноза экологического состояния заданного региона.
1. Уравнение нестационарного переноса загрязняющих примесей, распространяющихся в приземном слое атмосферы в условиях турбулентности, с начальными и граничными условиями имеет вид [1]:
А А А
д а( t Х") А А А А Я А А А А А А
’ -+a(t )-q(t ,х)+ р— (V(t ,x)-q(t ,х))~
д
~Гto
д t дх
А
8 q(t ,х)
К(t,х)-
д х
(1)
аА ААААа АААа А А _ _ А г ,
q(t = 0, х) = q 0(х), ,х = §) = ), q(t, х = I) = q■^(t), Г е [0,1] х е [0,1]
Все функции и переменные в уравнении (1) нормированы. Исходными данными уравнения (1) являются функции скорости ветра ¥(! , 'х), турбулентной диффузии К(1 ,х). источника примесей 8(1, х). коэффициента а(1). значения которых определены путем измерений, хранятся и обрабатываются в информационно-измерительных системах. Требуется построить вычислительную модель уравнения (1) для получения значений неизвестной функции
ц(1 ,х) - концентрации загрязняющих примесей в приземном слое атмосферы.
Построение вычислительной модели включает в себя следующие этапы. Для аппроксимации функций q(t,x), I ’(1,х). К(1,х). 8(1,х). представленных дискретными измерениями, а также их производных, применяется метод конеч-
ных элементов [2] с соответствующим выбором базисных функций. Выполняется редукция (1) к системе линейных дифференциальных уравнений первого порядка на основе метода взвешенной невязки относительно коэффициентов у\.(1).(\.(1)}. Затем система линейных дифференциальных уравнений первого порядка приводится к системе линейных алгебраических уравнений по схеме Кранка-Нико леона [4]. Неизвестные коэффициенты {Г )■(!) |. к = О, п +1 определяются путем минимизации функции невязки соответствующими методами [5]. Ниже подробно излагаются этапы построения вычислительной модели, обсуждаются вычислительные аспекты методов и алгоритмов.
2. Выбор и свойства базисных функций в методе взвешенной невязки. В соответствии с методами взвешенной невязки и конечных элементов искомая
функция д^,х) = ц(1 , х) представляется полиномом конечной степени (/„(1,х):
д„ (/, х) = С0 (/) • и0 (х) + 2 Ск (0 • ик (х) + Си+1 (0 • ми+1 (х),
к=1
где {ик(х)}„ = {///. (х)]и - совокупность базисных функций, определенная на
отрезке хе|Хо.Л | = |0.11. х = х; (С*(/)}, к = 0,и + 1 - коэффициенты представления функции д(1.х) в базисе {ик(х)}„ в момент времени I = I е [ОД].
В [1] для базисной функции {ик{х)}п = {м*(х)}и предложено следующее аналитическое представление:
Xх)=
Ґ к \а* ( к
Чк
х-х1
хк -хк \х2 х\ )
х9 -х
\хх х\ у
(2)
к к ”+1 где ик(х) определена на конечных элементах І1к = [х1 ,х2 ] \/хеО,к, ^ = и £\ ,
к=0
О = [0,1]. В этом вьфажении коэффициент дк характеризует амплитуду базисной функции ик(х); («/, Д) - некоторая система чисел таких, что ак > 0, Д > 0, к = О.я +1; Ик = | Х|/і. х* | - носитель к-й составляющей в аппроксимационной форме
п+1
и П+і(С(і),х)= ^ Ск^)ик(х)
Для аппроксимационной формы 1’„ , (('(О-'О совокупность {///(.г) | играет
П+1
роль базиса. Функция ( \(С(1).х) определена на множестве О = , О.
к-О
[ОД].
Базисная функция ик(х), соответствующая представлению (2), обладает рядом свойств. Основное - ее положительность Ух е | х*. .V2 ], т.е. ///.(.г) > 0. Другое - унимодальность, т.е. каждая из функций базиса {ик(х)} имеет единственную точку максимума на своем носителе 0.к.
* ^кх2 Ркх1 /04
хк = 7 , (3)
ак+рк
где х* < х*к < х% . Положение х*к на С1к определяется значениями чисел (ак. Д). Соотношение (3) определено для всех к, и поэтому для базиса {ик(х)} может быть построена последовательность {х|: . __г* |. Считая, что все пары
чисел (ак, Д) различны, приходим к тому, что среди точек х*к нет совпадаю-
щих. Совокупности |с*,х2,...,х*], 0.к- исходные характеристики базиса. Кроме того, поскольку каждая составляющая функция базиса {ик(х)} характе-
%
ризуется, помимо носителя О*, точкой ее максимума хк , пара чисел (ак, Д), связанная с ик(х), уже не может быть произвольной и ее выбор подчиняется условию
Рк=ак
( к * Д Х2 ~Хк
* ІГ
Кхк~х1 ,
(4)
Соотношение (4) показьшает, что параметры базиса (ак, Д) - вещественные числа, но могут быть и целыми.
Возможны различные варианты выбора характеристик базиса
[г*к,£1к,(ак,Рк)\. Если ввести вектор £ = {ь\}, к = 0,и + 1, где як =хк и <
Л/. < Л/, і, то получим соответствующий базис
Очевидно, что выбор 5= } однозначно определяет покрытие 0.к и («/.,
Д), т.е. определяет свойства базиса в целом. Подобный принцип построения базиса {ик(х) | приводит нас к аппроксимационному аналогу I 'п+\ (х. С. 3). определяемому не только вектором коэффициентов разложения С, но и вектором параметров £
Последним параметром исходной аналитической формулы (4) является величина дк, выбор которой определяется способом нормировки базисных функций ик (хк) = 1 или /.//.(.V/..) = 1:
Обозначим А\= sk- sk , и As2= л,, , - sk. В общем случае А\ Ф /£2. Эго
означает, что расположение узлов {sk} на интервале [0,1] может быть произвольным. Желательно, чтобы выбор S, определяющий характеристики базиса {sk, Q.k, (ак, Д)}, ставился в зависимость от общей структуры аппроксимируемой функции fix), хе[0,1], аналитическая сложность которой характеризуется ее вариацией. В [1] показано, что локальное поведение Un+1(x,C,3) может отражать особенности локального поведения fix) для xeQ.k при соответствующем выборе S и (ак, Д). Для этого разработаны алгоритмы «оптимизации» параметров базисной функции с учетом структурной сложности аппроксимируемых функций. Алгоритмы позволяют соответствующим образом располагать узлы на интервале [0,1] и для каждого С1к выбирать пару чисел (ак, Д). При этом, если Ф А'\ , то ак Ф Д. Тогда базисная функция ик(х) на носителе Q.k будет несимметричной. Расчеты показали, что применение алгоритмов «оптимизации» параметров базисной функции в аппроксимации fix), хе[0,1] с учетом ее вариации позволяет получать результаты с приемлемой точностью, не увеличивая при этом количество узлов. Это особенно важно для сложных функций.
Для равномерного распределения узлов {л/ ! на интервале [0,1], когда
= А2 = А, имеем следующий вариант аппроксимационной формы Un+1(x, С, 3), получаемый на основе (2), (5) и (6):
ик(х) =
-(51-х)а, к = О
------(х - яп )“ , к = п + 1
А“
где хе£1к= [^_ъ ^+1], О. = и£4 , = [0,1], ак = Д = а.
£=0
Поскольку ак = Д = а, то базисные функции %(х), хеИк будут симметричными. Кроме того, параметр а можно варьировать. Его значение влияет на «ширину» базисной функции. Чем больше значение а, тем «уже» йк (х). и наоборот, чем меньше значение а, тем «шире» йк(х). Варьировать значение а целесообразно с учетом вариации аппроксимируемой функции.
Третьим способом задания базиса является случай, когда 0.к = [0,1]
У к = 0, п +1. Тогда приходим к аппроксимационной форме II п+^ (х, С, ■§):
ик(х):
(1-х)
п-к
к = 0
- к гл \П-к
9кх 0-*) хк, к = п +1
£ = 1,п •.
(7)
где хеО* = [^_ь *м] = [0,1], £1= и£4 , = [0,1],
к-0
^ Д2 , акФ Рк, ак = к, Рк = п - к,
Чк =
к = 0 1
1, к = п + 1
к = 1,
У£ = 0,и + 1.
Если положить в (7) дк = , то получим (х) = рпк (х), £/„+1 (х, С, ^ =
5„+1(х)- многочлены Бернштейна. Таким образом, исходная аппроксимаци-онная форма £/и+1(С(/),х) является частным случаем многочлена Бернштейна. Это позволяет утверждать, что 1(С(1).х) наследует все свойства многочлена Бернштейна [3]. При этом она обладает дополнительными преимуществами, поскольку позволяет варьировать параметры базиса с учетом структурной сложности аппроксимируемой функции и, таким образом, повышать точность аппроксимации без значительного увеличения количества узлов.
Поскольку вычислительную модель исходного уравнения (1) предполагается строить на равномерной сетке узлов {х,}, / = 0,т +1, то {.V,! = = {х,} и для аппроксимации исходных данных в уравнении переноса нужно использовать
форму £/и+1(х, С, Б) с базисной функцией ик (х) = ~йк (х). Если же вычислительная модель строится на неравномерной сетке узлов, то в качестве базиса необходимо выбирать форму 11п+1(х, С, $ с базисной функцией ик(х) = ик(х). Для функций I '(/,х). /\(1,х) и 8(1,х). представленных дискретными значениями и входящих в уравнение (1), аппроксимационные выражения получат следующее представление:
Щ,х) * Го (I/)+ 2^^ (1} X (х) + Vт+х (/;) = Ут+1 ,х), (8)
2=1
Щ , *) * Ко «}) + Кг (/; )щ (х) +кт+1 (/;) = Кт+1 (/; , X) , (9)
2=1
£(/; , х) и ^0 (/; ) + (/; )мг (х) + £т+1 (/;) = 5'т+1 (/;, х) , (10)
2=1
где ¥г «} ) = Щ , X, ) , К, (I} ) = Щ , X,- ) , 5; (/; ) = £(/; , X,- ) .
3. Построение вычислительной модели на основе метода взвешенной невязки. Для исходного уравнения переноса (1) применим метод взвешенной невязки [2] и получим уравнение:
|ю,(х)
8 ф , X) + А («) Л (А^ А ) + . А ^ ^
5х 5х
Л
Г(Г,Х>
<Эд(Ч,х)
б/х = 0.
в
котором (®г(х)}, / = О./м +1 - система весовых функций, хе ^l/: О/ - носитель
т+1
®/(х), = О. О = [0,1]. Это выражение равносильно следующему:
г=о
5 q(t, х)
5/
А + «(/) ja>l(x)q(t, х)ск +
их
с!х -
(11)
дх дх
сЬ - £, |(х)5 (/, х)ск = 0,
Л Л
в котором обозначено / = / , х = х. Каждое слагаемое в выражении (11) преобразуем, используя аппроксимационные формы:
q(t,x) = ^ck(t)uk(x),
к-О
т+1
^ Ш, х)) = уск (t)uk (х), dt to
т+1
-ггШ’х)) = У! ск УЖ (*),
ох
(12)
(13)
к-О т+1
Й(^ х)) = 'Уск (/)< (х),
/и+1
где ик(х) = йк(х), хеО/. = |л> .V/. 11. О = и^Л; , ^ = [ОД]. Выполним преобра-
к=0
зования последовательно для каждого слагаемого (11).
Первое слагаемое (11) с учетом (12) получит вид:
J 0>i(x)
5 ^(f, х) н S II m+1 . I С,*(Г)мЛ(л)
dt Q L = 0
dx =С о (?) • J а>i (х)и о {x)dx fii.o
+ Е С^(/) | ®г(х)м^(х)А: + Ст+1(/) |®/(х)мт+1(х)А:,
^ = 1 ^,т+1
где О/ /- = 0.1 Г;(}к. Обозначим
п1,к = \х?1 , ] >1 >к = 0’т+1 > = I ■т1 (*)(*>* •
Тогда окончательно получим
А
д </(Г,х)
J ®;(*)
& =аг>0С0(Г) + T.alkCk(t) + al
к=\
(14)
(15)
Второе слагаемое (11) преобразуется аналогично первому с помощью соотношений (13) - (15). Имеем
cc(t) j ml (х) q(t, x)dx = a(t)(al 0C0 (t) + £ al kCk (t) +
^l,m+1 ^m+1 (0).
Q, fc=l
Выполним преобразования третьего слагаемого (11).
Р- I ®/(*)
п.
—(V(t,x)q(t,x))
ОХ
сЬ =
(17)
= Р ■ | соl(x)V'(t,x)q(t,x)dx + Р ■ la>l(x)V(t,x)q'(t,x)dx. п, п,
Интеграл в первом слагаемом выражения (17) можно преобразовать:
Л Л Л
| а>і (х) х) ції, х)ёх = Со (ґ) | (х) У(і, х)м0 (х)й?х +
0.1 &і0
т л л
+ 2 С* (0 1 ®/ (х) У'(1:, х)м^ (х)й?х + Ст+і (/) |®/(х)К'(/,х)мт+1(х)^х. (18)
^ = 1 &1,к &1,т+1
Для каждого интеграла, входящего в вьфажение (18), выполним процедуру интегрирования по частям. Окончательно получаем вьфажение для третьего слагаемого:
Р ■ I ®/(*)
Л = Р ■ (С0 (/)й, о (?) +
+ 2 С* (/)6,* (?) + Сга+1 (/)*/>т+1 (/)),
Ь/,о(0 =
У(1,Хд2д)~ \У(1,х)а>'[(х)и0(х)с1х, 1 = 0, к = О
п,„
- $У(1,х)ю'[(х)и0(х)сЬ, 1 = \,т+\, к = 0,
П,.о
(19)
Ъ1к({) = - \У{1,х)ю’1{х)ик{х)ёх, 1 = 0,т+1, к = \,т ,
^І.І
\У{ї, х)т'1(х)и т+1 (х)с1х, 1 = 0, т, к=т+1
Ь1,ш+1 (О =
■1/(Г=Хт+1,т+1)“ |1/(Г=*)®К*)«т+1 (*)*, 1=ГП+\, к = ГП +1 .
^1,т+1
Применение интегрирования по частям в выражении (18) позволило изба-
Л
виться от производной функции скорости ветра У'(1,х). Теперь в выражениях
Л
для коэффициентов {ЬцХ^} функцию I (1.x) можно заменить её аппроксима-ционным аналогом ¥т+\(1.,х) (18).
Выполним аналогичные преобразования четвертого слагаемого выражения (11). Имеем:
д х) дх
(Зх =
= у • \ \ о) 1(х)К'(t,x)q'(t,x)dx + \ аз 1(х) К ^ ,х) д" ,х)<&\.
[О; О/
Преобразуем первый интеграл данного выражения.
| а1{х)К'{1,х)^{1,х)ёх = | ®г(х)АГ'(^х)
о, о,
т+1
ТСк(!)и'к(х)
к-О
ёх =
= С0(0 | со1(х)К'(1,х)и'{)(х)с1х+ 2С^(0 \со1(х)К'(/, х)и^ (х)й?х+
О/ п к=\ О,
+ Ст+1(0 1а>1(х)к'(^х)и'т+лх)с1х-
^1,т+1
Для каждого интеграла первого, второго и третьего слагаемых выражения (11) выполним процедуру интегрирования по частям, после чего получим соотношение для четвертого слагаемого:
г- I ®гО) и,
КЦ,х)
8 х) 8х
(Ьс =
Г ' 1 со (0^1, о (О + ЕС ^ (0<*1к (О + С„ + 1 (Г )<*Л„ + 1 (О
^/,0 (0 -
-----К((,Хдд)- \К(1,х)со'1(х)и'0(х)сЬс, 1 = 0, к = 0,
А П,п
- $К(1,х)ю'1(х)и'0(х)сЬ, 1 = \,т+\, к = 0,
П
(20)
с11к(1) = - \К{1ух)со'1 (х)и’к(х)сЬс, / = 0,т+1, к = 1,т
о,.
^1т +1 (О
- (/, д:)й?/(х)м^+1 , 1 = 0,т, к = т + 1,
- — т+1) - / = т + 1, к = т + 1.
д ’ п,„,
В выражениях для коэффициентов {с//,/.(0 | также отсутствует производная
Л
для функции К (1.x). Для неё теперь будем использовать аппроксимацион-ную форму Кт ](1,х). определяемую вьфажением (9).
В пятом слагаемом (11) обозначим
gl{t)= \o)l{x)S{t,x)dx (21)
и для функции источника 8(1. х) будем использовать аппроксимационную форму 8т1(1.х). определяемую вьфажением (10).
Завершим преобразование (11), подставив в него (15), (16), (19)- (21). В результате получим:
т
а1, оСо(?) + 2 а1,кС к(г) + а1,т + \С т + \(1) +
к = 1
+ а(/)|а;>0С0(/) + ^а1кСк(1) + а;>т + 1Ст + 1(/)| +
+ Р-\с0 (Щ0 (0 + Д С, (Щк(1) + Ст+1 (Щт+1 (о| - (22)
- т • |0) т 0 (о+д ск а)(11к (?)+ст+1 (о^,т+1 (о| - # • & (о=о.
Систему (22) линейных дифференциальных уравнений относительно неизвестных коэффициентов К '/ (0!. к = \,т можно преобразовать к виду
т . т
2 а1кСк (0 + 2 ги (ОС* (0 = к, (!), (23)
к-1 £=1
если обозначить ги(/) = «(О^,* + Р■ Ъ1к(0~у-й1к(/),
Ь (0 = # • Яг (0 - «г,оО) (0 “ ^,^+1^+1 (0 “ гг,о (0' О) (0 - и+1 (/) • Си+1 (/) =
Л Л
= # • gl (0 - «/,(А (0 - «/,^+1^+1 (0 - г/;0 (0 • (0 - г/т+1 (0 • дт+1 (0,
1,к = 0, да + 1.
Систему (23) можно записать в матричной форме следующим образом:
А С = -Л С + к , (24)
где А = {ацс}, Я = (г/Д/)}, С = (С*(/)}, А = {А/(0}, 1,к = 0,да + 1. Как видно из (14), выражения для определения элементов матрицы /4 полностью определяются значениями весовых и базисных функций (®г(х)} и {ик(х)},
1,к = 0, да + 1 и не зависят от временного параметра I. В то же время, матрица
Я и вектор к должны вычисляться для каждого фиксированного момента времени. Кроме того, очевидно, что вектор к определяет поведение системы (24) на границах интервала [0,1] и включает в себя граничные условия, задаваемые для исходного уравнения (1).
Для модели (24) важным является выбор системы весовых функций (®/(х)}. Существует несколько подходов к их выбору [2]. Например, можно привлечь дополнительную физическую информацию. Если это невозможно, то можно положить ®/(х) = щ(х). При этом мы придем к системе дифференциальных уравнений, решаемых по методу Галеркина [2].
Другим важным моментом, возникающим при решении системы (24), является вычисление пределов интегрирования области = О/ г£ 1к. которые определяются согласно алгоритму, подробно изложенному в [1].
4. Вычислительный метод решения системы А-С = -И-С +к . Решение системы линейных дифференциальных уравнений (24) основано на построении вычислительных алгоритмов, приводящих её к системе алгебраических уравнений. Численное решение системы (24) осуществляется на основе ее
дискретизации по параметру / = / е [од]. Для этого можно использовать метод Кранка-Николсона [1, 4], являющийся одним из методов построения системы разностных уравнений для (24). Согласно этому методу, система (24) может быть записана следующим образом:
С (1 ] + 1 ) + С(1])
>Уг, (25)
}+
у2 _ +Я3
Г = /у+1 — 7 ■, /,к = 0,т +1, ] = 0,п -1.
Система (25) содержит неизвестные коэффициенты С1+1 = {с*+1 ], к = \,т , для определения которых выполняется построение функции невязки:
„Г - - »Г'А-
с; =9»° =<7(»о =(>.*»); ' Ч: С^1=й»1 =
= q(tj,xm+l); к = \,т , I = 0,т + \, ] = 0,п-1.
находятся путем минимизации функции невязки
. . 1/ , 1 т +11 . , |
Ш(С^)= г_!_.Т(рГ)2 или^(С^) = —-• I \рГ1\
\т + \ ыо 1 ) т +1 1=0 I I
Окончательно решение исходного уравнения (1) имеет вид:
чИ\ (*) = 4о;+1 + I СГ1 (0 • ик (х) + дС+\ , У = О,п -1, Ухе [ОД].
к=1
5. Заключительные замечания. Особенностью данной вычислительной модели является использование в ней аппроксимационной формы ( ](С(1).х).
Это означает, что для каждого фиксированного момента времени функции <7(/, х) ставится в соответствие вектор коэффициентов С(1) =
= {Ск(!)}, к = 1./м , являющийся представлением решения в базисе {ик(х) | в узловых точках хе[0,1], к = О./м +1, Определение значений неизвестных коэффициентов С(1) обеспечивается методами минимизации функции невязки. С помощью аппроксимационной формы ит+\ (С(/),х) получаем непрерьшный аналог искомой функции (/(I, х) \/хе[0,1] в каждый фиксированный момент времени I е[0,1].
В вычислительной модели используется новое аналитическое представление базисной функции. На его основе предлагаются различные способы формирования базиса и соответствующих аппроксимационных форм, таких как ГИ( |(х.С(1).,Ч) для случая неравномерного распределения узлов хк на всем интервале [ОД]; и т+1(х,С(!),8) для равномерного, 0т+1 (х,С,£) для случая, когда ит+1(х,С(!)) = Вт+1(х,С(!)). Таким образом, распределение узлов хк на интервале [ОД] может быть неравномерным. Кроме того, аппроксимационная форма (х.С'(О.Л') при определенных условиях может стать многочленами Бернштейна. При этом наследуются все свойства этих многочленов. С другой стороны, можно варьировать характеристики базиса {як, £ 1к. (ак, Д)}, к = О.т +1 с тем, чтобы учитывать структурную сложность аппроксимируемой функции. Здесь необходимо отметить, что базисные функции {ик(х)} принимают значения на конечных элементах \/хеО,к, к = О./м +1. Все эти дополнительные свойства и возможности аппроксимационной формы , (х, С.8)
позволяют аппроксимировать исходную функцию по ее дискретным отсчетам значительно лучше, чем многочлены Бернштейна, и при этом требуется меньшее число узлов. Однако есть и недостаток. Он состоит в том, что в отличие от многочленов Бернштейна, представление Ст \ (х,СЛТ) Ф ф/’(х). т.е. нет возможности аппроксимировать производную функции. Для решения этой проблемы можно привлечь технику конечно-разностной аппроксимации производных, однако это в свою очередь потребует значительного увеличения количества узлов и, как следствие, приведет к накоплению вычислительных ошибок. Благодаря использованию при построении вычислительной модели техники интегрирования по частям удается уйти от необходимости аппроксимировать производные функций исходных данных. Поскольку построение
вычислительной модели основано на методе взвешенной невязки, то требуется осуществлять набор весовых функций, к выбору которых существуют различные подходы. В нашем случае в качестве весовых функций выбираются сами базисные функции, что соответствует методу Галеркина. Но есть возможность использовать и другие весовые функции и, таким образом, расширять возможности вычислительной модели. Вариационный подход реализуется путем применения методов минимизации функции невязки, что позволяет уйти от непосредственного решения системы линейных алгебраических уравнений и связанных с этим известных проблем. Здесь необходимо заметить, что в вычислительном алгоритме модели применяются методы минимизации функции многих переменных нулевого порядка [5], позволяющие минимизировать функцию модуля невязки. Это является еще одним положительным свойством предложенной в данной работе вычислительной модели нестационарного уравнения переноса.
Литература
1. Семенчин Е.А., Наац В.И., Наац Н.Э. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.
2. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М., 1986.
3. Микеладзе Ш.Е. Численные методы математического анализа. М., 1953.
4. Мак-Кракен Д., Дорн У. Численные методы и программирование на фортране. М., 1977.
5. Дэннис Дж., мл., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М., 1988.
Северо-Кавказский государственный технический университет 16 апреля 2004 г.
