CC BY
13
УДК 519.5
ОПРЕДЕЛЕНИЕ ПРОИЗВОДНЫХ ЭМПИРИЧЕСКИХ ФУНКЦИЙ МЕТОДОМ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ В ЗАДАЧАХ ПЕРЕНОСА
© 2005 г В.И. Наац
The calculation methods for estimating derived empirical function, which is using for task solution of equation of the pollutants transfer in the frontier layer of the atmosphere, is described in this paper. The methods of numerical differentiation based on reduction this operation to the solution of first-genus Fredgolm's integral equation.
Эффективное применение численных методов теории переноса в задачах оперативного прогноза распространения загрязнений в пограничном слое атмосферы требует сопутствующей оперативной оценки характеристик атмосферы. Последнее обстоятельство требует в свою очередь разработки численных методов обработки и интерпретации экспериментальных данных, а также решения сопутствующих обратных задач теории переноса субстанции в турбулентных средах. Известно, что коэффициенты турбулентного обмена в направлении координатных осей в условиях пограничного слоя атмосферы могут быть функционально связаны с производными компонент векторного поля скорости ветра по пространственным переменным. В статье рассматривается это направление исследований, сформулированы соответствующие модели, для которых исходными дан-
\дУ1 —
ными являются матрицы <—- >, i,j = 1,3, где Vi - компоненты вектора
[dxj
скорости ветра V(P,t) в приземном слое атмосферы. Численное решение подобных математических задач требует алгоритмов дифференцирования эмпирических функций, построение и исследование которых далее подробно излагаются в статье.
1. Вывод интегрального уравнения для численного дифференцирования эмпирических данных. Полагаем, что функциональная зависимость переменных x и y (тоже y=f(x)) представляется совокупностью {(xk,yk )}n, в которой значения yk измерены (известны) с относительной ошибкой Sk, не превышающей некоторого значения s (sk < s для
Vk = 1,n). Подобные функциональные зависимости можно условно называть эмпирическими функциями. При решении уравнений переноса [1] в качестве подобных функций выступают исходные данные, такие как поле скорости ветра и поле турбулентности, определяемое коэффициентами турбулентной диффузии, и представленные измерительными данными {v(Xk,tj )} и Xk,tj )} соответственно. Методика оценки производных эмпирических функций представляет собой задачу, которая относит-
ся к так называемым некорректным задачам математического анализа и требует соответствующего подхода и соответствующих численных методов. Далее излагается метод численного дифференцирования, основанный на сведении этой операции к решению интегрального уравнения Фред-гольма первого рода.
Пусть задана функция /(х) в пределах интервала [х0 ,хп ]. Требуется
вычислить /'(х). Поскольку связь /(х) с её производной /'(х) представлена через первообразную, то может быть записано следующее исходное выражение
/(х) = /(Хо) +5 /или 5 /'(Г)Л = /(х) - /(Хо). (1)
х0 х0
Введем вспомогательную функцию
г1, если х0 < / < х К(х, ) = \ 0 (2)
[0, если х < / < хп
и перепишем (1) в виде
} К (х, Г) /' (/)Л = / (х) - / (хо). (3)
хо
Формально (3) является интегральным уравнением относительно неизвестной функции /'(х) , а правая его часть определяется исходной функцией /(х). Заметим, однако, что в соответствии с (2) ядро К(х,) не является всюду непрерывной функцией в области своего определения £2 = \хо,хп]х \хо,хп], поскольку оно имеет скачок в точке / = х . Это обстоятельство затрудняет практическое использование (3), так как в соответствии с теорией интегральных уравнений гарантировать существование и единственность решения можно только в случае непрерывных ядер. В связи с этим необходимо (3) трансформировать в интегральное уравнение с непрерывным ядром, которое при этом оставалось бы эквивалентным исходному. Подобную операцию можно выполнить, прибегая к оператору К*, сопряженному к исходному интегральному оператору К, связь между которыми определяется следующим соотношением:
(Ки, V ) = (и, К *у ), (4)
где и и V - некоторые интегрируемые функции по переменной х.
Найдем ядро К *( х, /) сопряженного оператора К*, исходя из (2) и (4). Имеем для левой части (4):
J J K(х,t)u(t)dt v(x)dx =J u(t) J K(x, t)v(x)dx
dt. (5)
Учитывая (2) и сопоставляя (4) с (5), устанавливаем, что в качестве ядра сопряженного оператора следует принять функцию * 10, если х0 < х < *
К * (х, *) = \ * < < . (6)
[1, если * < х < хп
С учетом соотношений между переменными х и * в (6) ядро К *( х, *)
*
следует писать в виде К (*, х) по аналогии с записью исходного ядра
К(х,) в (3). Далее удобно ввести вместо I новую переменную 5. В
итоге ядро сопряженного оператора примет вид
* го, если х0 < х < 5
К *(5,х) = \' 0 . (7)
[7, если 5 < х < хп
Оно разрывное на диагонали в точках 5 = х, как и исходное К(х,). Построим новую функцию двух переменных К (5, *), положив
хп
К(5,0 = | К (5, *)К (х, *)йх. (8)
х0
Вычислим интеграл (8) для двух ситуаций:
а) х00 < * < 5 < хп;
б) х0 < 5 < * < хп .
Для случая а) интеграл в правой части (8) может быть представлен в
виде суммы трех интегралов, а именно
хп » * „
I К (5, *)К(х,*)йх = | К (5, *)К(х,*)йх +
х° хо (9)
5 хп
+1К*(5, *)К(х, *)Сх + I К*(5, *)К(х, *)Сх.
* 5
С учетом (7) и (2) выражение (9) получит значение (хп - 5). Для случая б) будем иметь
хп „ 5
I К (5, *)К(х,*)йх = I К (5, *)К(х,*)Сх +
х0 х0 * хп
+1К*(5, *)К(х, *)Сх + I К*(5, *)К(х, *)йх.
5 *
Прямое вычисление этого интеграла дает значение (хп - *). В итоге получаем
гхп - 5, если х0 < * < 5 К (5, *) = \ п 0 . (10)
[хп - *, если 5 < * < 5
Непрерывность этой функции на диагонали 5 = * очевидна. Остается данное ядро ввести в рассматриваемую задачу (3). С этой целью умножим
(3) на К (s, х) и проинтегрируем в пределах от х0 до хп по переменной х. Имеем
х {х х х
|К*(?,х) | К(х,г)dx = | К*(5, х)/(х^х- /(х0) | К(5, х)ск . (11) х0 v х0 / х0 х0
Меняя порядок интегрирования слева в (11) и учитывая (8), приходим
хп ~ хп *
к интегралу | К(5,г)/'(г)dt. Справа в (11) | К (5, x)dx = хп - 5, и, значит, х0 х0 имеем вместо (11) выражение
хп
1к (5, г) / '(^ = ^), (12)
х0
хп
в котором обозначено = | /(x)dx - /(х0)(хп - 5).
Полученное интегральное уравнение (12) относительно /' (г), г <е\хо ,хп ] является интегральным уравнением Фредгольма первого рода с непрерывным ядром. Его численное решение можно построить на основе вариационной задачи для так называемого сглаживающего функционала, что ведет к регуляризирующему оператору К^1. Численно он реали-
- 1 I Т I 1 т
зуется с помощью матрицы Ка = 1К К + аП К , где К - матричный
аналог исходного интегрального оператора К в (12). В итоге искомое решение определится так:
Га=(КТК +а1 У Ктф . (13)
2. Вычислительная схема метода, результаты её численной реализации. Предположим, что распределения У(х,г) и V'(х,г) заданы в «Блоке исходных данных» функциями
П | | (, „ „„ ■ (П
V(x,t) = V0 1 + 0,25sinj — xjj^1 + °,85•sint
V'(x, t) = V° • 0,25 • — • j cos x j j • + 0,85 • sin t
значения которых вычисляются в точках равномерной сетки (xj-,tj) и представлены массивами {V(xi, tj )} и {vt (xi, tj)} = {V'(xi, tj)}, i = 0,m,
j = 0,n . Заданы массивы {xi}, i = 0,m, Vxi e [0,X]; {tj}, j = 0,n , ytj e [0,T]. В этом случае дискретизация задачи (10), (12), (13) примет вид
I xm - xi, если x0 < xk < xi - -
К(xi, xk) = < , k = 0,m , i = 0,m ;
xm - xk , если xi < xk < xm
f(xk) = V(xk,tj), f'(xk) = VT(xk,tj), к = Ö~m ; j = Ö~n ;
m
тщК (x;, xk) f'(xk )ax = p( x;), (14)
k=0
m -
ф(xi ) = z ®kf (xk )Ax - f (x0 )(xm - xi x i = 0,m ,
k=0
a>0 = com = 0.5 ; щ = 1, k = 1,m -1. Выполним нормирование переменных и массивов:
xi ■ -— Л tj
xi =—, i = 0,m; t: =—, j = 0,n ; ' x> i t
V (xkt:) л л л V' (xkt:) -
V (xk, tj) = —, V •( xk, tj) =—kj , k = 0,m , j = 0,n . V0 V0 Переходя к нормированным величинам в выражениях (14), получим
К(xi, xk) =
СС (xi, xk) = X • К (xi, xk), 1 - xi, если 0 < xk < xi
k = 0,m , i = 0,m,
1 - xk , если xi < xk < 1
X • Vo • z щ К (xi, xk) f'(xk) a x = X • Vo •
k=0
i = 0,m,
или
z®kf (xk )a x- f (0)(1 - xi)
k=0
Л АЛ АЛЛ
Z Щ К (xi, xk) f (xk) a x = cp( xi),
k=0
<р(хг) = ^ ак /(хк )д х- /(0)(1 - хг), I = 0, т . (15)
к=0
Построение вычислительного алгоритма для определения значений <|/(хк)| основано на вычислении матриц:
1) А = ((КТК + а1), где К =а ■ К-Дх ;
2) А+=(ТК + а1) *(ЖТК + а1 )"7 - обобщенная обратная матрица, вычисляемая с помощью итерационного метода Бен-Израэля [2]. Тогда имеем приближенное решение
f a*(KTK + aI)+ KTv = f' .
(16)
Параметр регуляризации а определяется по невязке [3]. С учетом вышесказанного алгоритм получения приближенного решения (16) будет включать в себя следующие шаги: 1) ] = 0;
2) / (]) = Уг,}, I = от;
3) моделируется процесс появления случайных ошибок 8 при измере-
Л (8) я (8) я
нии «экспериментальных» значений р (]): р1 (]) = рi(]) • (1 + ©¿8), где ©г- - случайные числа, равномерно распределенные на интервале (-1, 1),
генерируемые
Л (S) Л
р (j) -р(j)
датчиком = °Ф(j), т-е-
случайных
чисел.
При
j) =
7-1) 2
(m +1) i=o
(Л (s) Л Рг (j) -Рг(j)
a<p=db) 2a<p(j )
(n +1) j=0
этом
(17)
где рг (]) =2 /к (]) ах - (/0 (])) • (1 - х), /к (]) = /(хк, ),
к=0
/0 (]) = /(х0, О), I = 0,т ;
4) задаются значения ао, 0 < г < 1;
5) р = 0;
6) вычисляется значение параметра регуляризации ар = аогр ;
7) вычисляется обобщенная обратная матрица:
а) С = К , С = {р,,к }, В = КТ , В = {р1Л }, Ь,,к = Ск, , к = 0,т, I = 0,т;
т ___
б) А = В х С, А = а к }, а1 к = 2Ь /С/ к , к = 0,т , I = 0, т ;
/=0
в) =
|ai , k + ap , если i = k
ai ,k.
если i Ф k
г) вычисляется A+ , D = A+, D = {di k }, k = 0, m , i = 0,m ; 8) приближенное решение f'a :
а) ° = \, 8г,к = Елг,1ъг,к, к = 0,т,г = 0,т;
I=0
т , №
б) Гар, (7) = Е Яг,к <Рг (7);
Р к=0
9) подставляя найденное решение в исходное уравнение (15), определяем величину погрешности
(
Ol(j) =
1
(m +1)
f m л _ л л № ^
fan ,k (j)Д *(j)
k=0 p
1
1_(n+1) jV1
j).
(18)
10) если ст1(]) <а1р(]), то а (]) = ар, перейти на шаг 11; иначе -р = р +1, перейти на шаг 6;
11) определяется значение погрешности приближенного решения
е( j) =
1
1
(n +1) j=0
Ss(j); (19)
12) % = а(j), г = 0,т;
13) ] = ] +1, если ] < и , то перейти на шаг 2;
иначе - конец алгоритма.
Численная реализация описанного выше алгоритма проводилась при следующих значениях исходных данных 7=600 с, Х=1600 м, У0=5 м-с-1. Диапазон изменения значений скорости ветра составил [Ут,„;Утсх] = = [5; 11,56] м-с-1, производной скорости ветра - [У'гш„;У'тса\ = = [1,56 Е-0,3; 3,63] м-с-1.
Первоначально выполнялись расчеты погрешностей с1 (18) при различных значениях а^ (17), зависящих, в свою очередь, от величины 5, с целью исследования влияния ошибок измерений в исходных данных с на вычислительный процесс. Так, если 5=0,01 (1%), то ст(!,=7,77Е-03, с1=3,57Е-03; если 5=0,05 (5%), то <тр=3,88Е-02, с1=1,75Е-02; если 5=0,1 (10 %), то а(р=7,77Е-02, о>=3,51Е-02. При этом в формулах (18) исполь-
г = 0,т,
зовались точные значения производной i / (j) i =1 V.
г. j
T \ Л
7 = 0,„ , генерируемые в «Блоке исходных данных».
Алгоритм решения обратной задачи далее реализуется при начальных значениях параметра регуляризации а0=0,25 и г=0,5. На рис. 1, 2 показа-
ны пространственно-временные распределения (профили) полей У'т и V' при значении 8= 0,01. При этом отклонение точного решения от приближенного (26) составило е=8,5Е-03. Профиль приближенного решения почти полностью повторяет профиль точного решения.
Проведено численное исследование влияния величины погрешности в исходных данных ар на точность получаемых приближенных решений. На рис. 3 показаны графики точного V'т и приближенных к нему решений V! , V', V' при значениях 8= {0,01; 0,05; 0,1} соответственно для фиксированного момента времени. В таблице приведены усредненные по времени значения величин е и аь а также коэффициент усиления ошибки получаемого приближенного решения от величины погрешности в исходных данных. Результаты расчетов показывают, что при наличии измерительных ошибок в исходных данных, составляющих до 10 % от их точных значений, вычислительный алгоритм решения данной обратной задачи является устойчивым (коэффициент усиления ошибки 0<к<1) и обеспечивает получение приближенных решений с приемлемой точностью е< 1,0Е-01. При этом для каждой строки таблицы выполняется соотношение а1<ар.
Исследование влияния величины погрешности в исходных данных и 8на точность получаемых решений еи а
5 Е k=gv!E
0 0 1,38Е-02 6,44Е-03 0
0,01 7,77Е-02 1,92Е-02 8,51Е-03 0,40
0,05 3,88Е-02 4,32Е-02 1,89Е-02 0,89
0,10 7,77Е-02 7,02Е-02 3,19Е-02 1,11
0,15 0,15 0,12 5,61Е-02 1,25
0,20 0,20 0,16 7,17Е-02 1,25
Рис. 3. Графики точного У'Т и приближенных решений V/, V', V3 при значениях
5={0,01; 0,05; 0,1} соответственно;]=6; е,(]=6)=1,89Е- 02, е2(]=6)=2,65Е- 02, ез(]=6)=3,22Е- 02
Проведенные численные исследования вычислительного метода определения производных эмпирических функций, основанного на решении интегрального уравнения Фредгольма первого рода методом обратной задачи, позволяют сделать вывод о вполне удовлетворительной, устойчивой работе соответствующего ему вычислительного алгоритма.
Литература
1. Семенчин Е.А. и др. Математическое моделирование нестационарного переноса примеси в пограничном слое атмосферы. М., 2003.
2. Кубанива М. и др. Математическая экономика на персональном компьютере: Пер. с яп. М., 1991.
3. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М., 1979. Северо-Кавказский государственный технический университет 20 апреля 2005 г.
