Расчетно-аналитические модели для уравнений параболического типа с приближенными данными на основе методов прикладного гармонического анализа и вариационного метода взвешенной невязки Текст научной статьи по специальности «Математика»

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

«НАУКА. ИННОВАЦИИ. ТЕХНОЛОГИИ». №3, 2015

УДК 517.972:519,633 Наац И. Э. [Naats I. Е.], Наац В. И. [Naats V. I.], Рыскаленко P. A. [Ryskalenko R. А.]

РАСЧЕТНО-АНАЛИТИЧЕСКИЕ МОДЕЛИ ДЛЯ УРАВНЕНИЙ ПАРАБОЛИЧЕСКОГО ТИПА С ПРИБЛИЖЕННЫМИ ДАННЫМИ НА ОСНОВЕ МЕТОДОВ ПРИКЛАДНОГО ГАРМОНИЧЕСКОГО АНАЛИЗА И ВАРИАЦИОННОГО МЕТОДА ВЗВЕШЕННОЙ НЕВЯЗКИ

Computational and analytical models for equations of parabolic type with approximate data on the basis of methods of applied harmonic analysis and variational method of weighted residuals

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

Ключевые слова: уравнение массопереноса, расчетно-аналитичес-кая модель, методы прикладного гармонического анализа, регуляризи-рующие численные методы, вариационный метод взвешенной невязки.

Considered non-stationary equation of mass transfer of impurities in the atmosphere, which is a linear differential equation of parabolic type. It is assumed that the mathematical model requires attracting approximately given initial data or semiempirical functions for their determination. To solve this problem in the work builds regulating numerical methods and related computational and analytical models based on the methods of applied harmonic analysis and variational method of weighted residuals.

Key words: the equation of mass transfer, computational and analytical model, methods of applied harmonic analysis, regularization numerical methods, variational method of weighted residuals.

При разработке численных методов и соответствующих расчетно-аналитических моделей для решения различных функциональных уравнений используются обычно те или иные методы теории прибли-

жения функций, в соответствии с которыми искомые функции могут быть представлены соответствующими обобщенными полиномами. Последующее построение решающего алгоритма сводится к нахождению более простых уравнений для вычисления коэффициентов аппроксимационных полиномов. Подобный подход в ряде случаев приводит к определенным затруднениям с точки зрения устойчивости вычислительных процессов, особенно тогда, когда исходные данные в задаче являются приближенно заданными функциями. Это замечание относится прежде всего к решению дифференциальных уравнений, поскольку операции дифференцирования на приближенных данных являются некорректными в математическом отношении. Аналогичные задачи возникают в математическом моделировании реальных явлений, описываемых дифференциальными уравнениями и требующих привлечения эмпирических данных или полуэмпирических расчетных формул. К подобным задачам можно, в частности, отнести задачу моделирования пространственно временной изменчивости поля концентрации загрязняющих веществ, распространяющихся в приземном слое атмосферы [1, 2, 3].

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

Постановка задачи.

Изложение предлагаемых ниже методов регуляризации осуществляется на примере уравнения параболического типа, каким и является так называемое полуэмпирическое уравнение переноса вида [1]:

й(х,у,г,1):

8_ дх

д

—(*»+-=-(*»+— СП»)

ду

д_ дг

+

+

дх ду дг

В этом уравнении (x,y,z,t) е ii = il + iî. u(x,y,z,t), - концентрация загрязняющих примесей в атмосфере, Vx(x,yzJ), Vv(x,y,z,t), K(x,y,z,t) - компоненты вектора скорости ветра V в направлении координатных осей Ох, Оу, Oz соответственно, K(x,v,z,t) - коэффициент турбулентной диффузии, S(x,y,z,t) - источник примесей, в точке пространства P(x,v,z) в момент времени t. При этом ù соответствует n't. Множество исходных данных в рассматриваемой модели в дальнейшем обозначим через В = {(V,K)} и считаем состоящим из функций, не дифференцируемых в обычном смысле. Подобное предположение естественно делает задачу решения (1) на множестве В неопределенной. Последующее изложение предлагаемого подхода к нахождению приближенных решений модели (1) на множестве В осуществляем в предположении, что эти решения описываются тригонометрическим полиномом

, . , ч 8 , ч . лпх . пту . жкг

u(x,y,z,t)^uNMK(x,y,z,t)= — —srn- —-sin----, (2)

abc n=1 m=1 k=i abc

определенным на брусс Q = iix x ûv x Q = [0, a] x [0, 6] x [0, с]. Полагаем далее, что u(x,y,z,t) на границе указанной области принимает нулевые значения. В противном случае требуется ввести соответствующие граничные элементы (составляющие компоненты в (2) [4, 5]). При этом функции базиса (в данном случае речь идет о разложении функций по произведениям), а именно

, ч 8 ппх . пту . 7vkz

qnmk (*» -V, Z) = —— sm - • sm - • sm--,

abc abc

в принципе могут быть и иными, например, это могут быть функции вида

8 ттх типу 1ÛZ

Рnrnk (*> У, Z) = ~Г C°S--008 к' ' COS-'

abc abc

а также их линейные комбинации. Если допустить, что функции V" (x,y,z, /) (/ = 1,2,3) и K(x,y,z,t) входящие в множество исходных данных В. представимы в виде полиномов (2), а именно VmIK (-v. у. z, t) и KNMK{x,y,z,t\ то возникает вопрос: в каком соотношении будут при этом

находиться их частные производные, скажем, V' и (Ум.ж)х' и т. д. в ситуации, когда данные являются приближенными? Для того чтобы ответить на этот вопрос, обратимся к некоторым известным положениям прикладного гармонического анализа, предварительно заметив, что понятие «прикладной гармонический анализ», используемое в настоящей работе, интерпретируется в соответствии с работой К. Ланцоша [6].

Положим для примера, что некоторая функция/(х) непрерывна на интервале [0, а] и представляется с приемлемой точностью конечным рядом Фурье по системе jsin (п = 1,2 ..), т.е. имеет место соотношение

ff /Itfx

/Ш а In О) = £а« sin — (3)

для всех х из указанного интервала. Подобное представление имеет практическую значимость, если \а„\ —> 0 при п —> да для данной функции /(х). Для непрерывной функции это действительно имеет место. Если же исследуемая функция/(х) представлена своим приближением /(х). содержащим помимо регулярной компоненты, для которой указанное условие выполняется, также и нерегулярные компоненты, для которых значение \а,\ не убывает с ростом п, то представление (3) теряет практический смысл. Напомним, что нерегулярные компоненты в исходных данных появляются всякий раз, когда используются в вычислениях так называемые эмпирические данные, содержащие ошибки измерений (реализации случайных процессов). Типичной ситуацией в задачах прикладного гармонического анализа является случай, когда \ a„(f)\ 0 при и —>■ да, и можно говорить о соответствии/(х) и/,(х) в смысле (3), однако эта сходимость последовательности j | а„ (/) J явно недостаточна для сходимости ряда ^a„(/)(sin . и, значит, нет соответствия между/'(х) и /'. V(.v). В связи с этим в дальнейшем будем различать последовательности коэффициентов Фурье \a„{f)} и {a„(f) = а„( /)}, которые при п —> / ведут себя различным образом. Заметим, что в каждой конкретной задаче смысл обозначения fix) и его соотношение с функцией/(х) определяются ее содержанием и требует четкой формулировки принимаемых предположений.

Задача суммирования продифференцированных рядов Фурье, коэффициенты которых вычислены по приближенным данным, и метод регуляризации. В контексте выше изложенного в прикладном анализе особое внимание уделяется задаче суммирования продифференцированных рядов Фурье, коэффициенты которых вычислены по приближенным данным. Приемы, связанные с регуляризацией поведения продифференцированных рядов Фурье, как правило, основываются на тех или иных способах их усреднения. Простейшим из них можно считать так называемый метод <7 -множителей, предложенный К. Ланцо-шем [6]. Не останавливаясь подробно на технике вычисления этих множителей, приведем итоговые расчетные выражения для производных от тригонометрических аппроксимационных полиномов типа (2). Для простоты ограничимся случаем функции двух переменных по х и у, а именно:

К О, У, 0 ~ иШх (х, апт(0 ■ с41 • [яп ппх

аЪс В=1 И=1 1

81П

пту

(4)

о n м

8 ТУя

ппх

й'у(х,у,0Яи'ш (х,у,^ = — (0•

аЬс^^ а

7(1) •

ят

пту

(5)

8Ш

,-0)

п п а' N

. ( п т л БШ----

и М

г п тл

Ъ М

(6)

Поскольку значения сг-множителей (6) удовлетворяют условию 0 < а < 1, то их роль в ускорении сходимости рядов (4) (5) вполне очевидна. Ясно, что при п = 0 = 1, и о^ = 0 при п = N - 1. Используя данный подход, нетрудно показать справедливость следующих приближенных выражений для оценки векторных производных функций и (лм'.г./). представленных своими приближениями / (тоже {а„ (у~) = а„ (/)}). Опуская промежуточные вычисления имеем:

п=1 т= 1

БШ-

ттх

•8111

пту

V

г

N М

л=1 т=1

n м

/

вш

лту

<т(2) =

БШ

я=1 т=1 ( К П

\2а ' ТУ

л

7ШХ

У

'я- ил

V V

2а ЛГ

ят

V «у

г(1)

'М,т

8111

типу

У У

Изложенный здесь метод регуляции тригонометрических аппроксимационных многочленов в случае применения к ним операции дифференцирования можно рассмотреть с более общей точки зрения, используя понятие оператора обобщенного дифференцирования, введенного ранее в работах авторов [2, 3]. Действительно, обозначим через &„(х) составляющую базиса {9„(х)} (п = 1,2). Тогда изложенный выше метод «дифференцирования» рядов приводит к соотношению:

(7)

где = —-ах

В соответствии с этим оператор /У следует считать оператором обобщенного дифференцирования для конечных рядов Фурье с приближенными коэффициентами ап(/) = ап( /'). Формально при N —>• / оператор 1)п> может быть сколь угодно близок к оператору I)1'' в

смысле подходящей нормы \\0 / /.) '/ || при условии, что ||/ - / || достаточно мала. Если речь заходит об обобщенных производных у-го порядка, то (7) можно писать в виде

= где -

/ 'п ча п 'лг, \

г \

п п

V и ■а N. У

тс

при 3„(х) = зш—пх . (8) а

Построение расчетно-аналитической модели для уравнения массопереноса. Выражения, подобные (8), могут быть построены не только для тригонометрических базисных функций. Используя представление искомой функции и(х,у, /) в виде тригонометрического полинома (2) и введенные выше правила его дифференцирования, можно вернуться к модели (1), определенной на множестве исходных данных В = {{у,К)} Предполагается, что приближенно заданные функции представлены последовательностями своих коэффициентов Фурье, т.е. для компонент вектора скорости ветра и коэффициента турбулен-

тной диффузии {апт (1,К)} соответственно. Для компактной записи соответствующих расчетных формул введем следующие обозначения

, ч 4 тх . пту

чЛх,У)=-г- sin--sin—-

ab a b ,

_4_ ab 4_ ab

, mx sin-

v a j

a

. типу -sin—— b

типу

sm

b

ab

4

/

же

(2,x)(

Чпт V ■*">?) , nm

ab

sin-

. mx sin-

v а у

7(1)

, шпу sm——

у ,

4

nnx

Фпш y\x,y) =—-sm

ab a

№ 'Mm

. типу

sin——

. nmy sm——

Исходную модель (1) удобно представить в виде ди ~ ~

— = -Q(x,y,t,u,V) + Q(x,y,t,u,K) + S(x,y,t\ (9)

dt

где первый функционал по и в правой части (9) соответству-

ет первому слагаемому в правой части (1), а второй - второму. С учетом принятых обозначений имеем

N М -

Q(x,y,t,u,V) = XZ2- (0 • Qnm(x,y,V\

п=1 »1=1

где

QJx>y>b=VU qnя(х,у)+ q™\х,y)+Vly-qJx,у)Щ^УК*,У1 V1(x,y)aV2(x,y)

- соответствующие компоненты вектора скорости ветра. Аналогично записываем и второй функционал в (9), а именно:

n м

Р(х, y,t,u,K) = YT."dnm (0 ■ рпт (*> У, К).

n—1 т—\

где

Рпт^уЛ)=К-д^{х,у) + K-q^\x,y) + К'/Ч^\х,у)+К ■ д?/\х,у) В итоге уравнение (9) примет вид

n м n м -.

££й(0 • q„fry) = -££гвт(0 • Bnm(x,y,V ,К) + S(x,y,t\ (10)

п=1 т-1 и=1 т=1

где Bnm(x,y,V,K) = Qnm(x,y,V)-Pnm(x,y,Ky

Уравнение (10) является линейной системой обыкновенных дифференциальных уравнений относительно матрицы коэффициентов {am(t)}, которые определяют Фурье-разложение функции и(х,у, /) в двухмерный ряд Фурье. Для их численного нахождения необходимо тем или иным способом «свернуть» систему (10) по пространственным переменным х и у. В частности эту операцию можно выполнить с помощью метода взвешенной невязки [1, 4, 5].

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

} wi\ ,i2 (*> У)р(х> У> t)dxdy = 0, (11)

п

гдеp(x,y,t) - невязка для системы (11) и /, = 1.....V. /2 = 1....А/.

Без ограничения общности можно полагать, что и (х,у) = шт (х) • у»и (у) и остается разумно выбрать весовые функции И' (.у), сообразуясь с теми или иными особенностями решаемой задачи. Соответствующие вопросы рассматривались ранее в работах [1, 4, 5]. В контексте данной статьи заметим следующее. Выбор системы весовых функций можно в определенной степени связать с аналитическими свойствами составляющих функций базиса. Использование в качестве последних тригонометрических полиномов Зп(х) = зт-^тс делает целесообразным по ряду причин в качестве Шу(х) также выбрать тригонометрические полиномы. Не останавливаясь подробно, уточним, что разумно в качестве Щ (х) выбирать многочлены Фейера [6] по тригонометрическим системам, т.е. положить

щ(х) = Ф

2лу

/

8111

ЛЛ2

V-л-

V

бш

л •

/ 7

где V - порядок соответствующего полинома. Выбор порядка V

согласуется с размерностью задачи N и величиной ожидаемой погрешности в исходных данных. С учетом того, что в интегралах (11)*^ Л (х, у) = ФУ1 (х, х1{) • Фу% (у, д ), где {(х^, у,2)} - соответствующая сетка в по пространственным переменным х и у, система (10) преобразуется к виду

А • а(0 = -В(0 ■ 5(0 + Л(0,

(12)

А =

ДфУ1 (х,хг1) • Ф„2 (у,у1г) ■ Чпт{х,у)(кйу

т=

Уравнение (12), представленное в виде =—А~]В(^)а^) + + А Н(р\ является уравнением эволюционного типа, методы решения которого ранее рассматривались авторами в работе [1]. Расчетные соотношения выше приводились для двухмерной задачи. Ясно, что не составляет труда расписать их для более сложного трехмерного случая.

В заключение, рассмотрим в качестве примера приложения изложенной здесь теории вариант одномерной задачи, акцентируя в большей мере внимание на способах регуляризации вычислительного процесса. Допустим, что требуется построить расчетно-аналитическую модель для нахождения приближенных решений уравнения вида

Л»

— + С0-и,х + С1-(и2)'х + С2-итххх=з(х, 0, т

(13)

в которой С о, С\ и С2 считаются постоянными. Последнее предположение для вышеизложенного подхода к решению уравнения (13) не носит принципиального характера, а лишь упрощает запись расчетных формул, приводимых ниже. Выражение (13) можно переписать в виде

й + р(и) ■ и'х + С2 • и* =

(14)

если ввести р(и) = (С„, + 2 С\ и). То, что в модели (13) присутствует производная третьего порядка, существенно затрудняет применение обычных разностных схем (сеточных моделей), особенно при наличии погрешностей в исходных данных (С0, С\. С2). Используя рассмотренные выше представления исходной функции и(х,1) и ее производных по переменной х е [0, а], записываем

. , ч 1 . , . . 7W.X ИN (JC,0= 2s ап (0 • Sin-

a S а

1 N а £i

sm-

*N,xxx\

a j,

sm-

a y

Значения о-множителей определяются в соответствии с (8). Уравнение (14) с учетом этих соотношений предстанет в виде системы уравнений относительно a{t)

£àB(0 ■ sin^ + р{а) ■ £a„(f) • a®,

f . mxл

sin-

а

\

-i-c2-£a„(o-<,

sm-

v a

-a-s{x,t)

(15)

Используя метод взвешенной невязки [1, 4, 5], систему (15) можно привести к стандартному виду

А ■ à(t) - -B(a,t) • a(t) + sit).

где элементы матрицы А определяются формулой

а

аы = fw,(x)-sin—¿х, I — \..N. n -I..N,

J n

(16)

bbi{t) = \wl{x)- p(a{t))'

о

матрицы />'(/) : f

, лт х sin—

а у

■a^ndx + C2]Wlix)-a%

f

. mx sm-

v

a

dx-

и вектора s (t) :

a

s jit) -a- w, (x) • s(x,t)dx

Применяя далее к (16) принцип рекурсий на сетке узлов {/ ■} по переменной t е [0,7], получим итоговое расчетное соотношение:

.дц± 0 5и _ 1) = _т)) аи , 1) - аи ~ О + = 1Д, (17)

2т 2т

в котором подразумевается, что а(/') = ), то же касается и ?(_/) Особенности вычислительных процессов, соответствующие расчетно-аналитической модели (17), подробно исследовались авторами в работе [1].

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Наац В.И., Наац И.Э. Математические модели и численные методы в задачах экологического мониторинга атмосферы: монография. М.: Физматлит, 2010. 328 с.

2. Рыскаленко Р.А., Чемеригина М.С. Операторы обобщенного дифференцирования в численных методах решения нелинейного уравнения переноса с приближенными данными // Вестник Северо-Кавказского федерального университета. 2013. №1 (34). С. 35-38.

3. Наац И.Э., Артемов C.B. Итерационные алгоритмы для численного решения уравнения переноса на основе операторов обобщенного дифференцирования // Вестник Северо-Кавказского федерального университета. 2013. №1 (34). С. 21-26.

4. Наац В.И. Вычислительная модель нестационарного уравнения переноса примеси на основе метода взвешенной невязки // Известия вузов. Северо-Кавказский регион. Естественные науки. 2004. №5(04). С. 3-15.

5. Наац В.И., Гаршина ТВ. Вычислительная модель нестационарного уравнения переноса примесей на основе метода взвешенной невязки и операторов обобщенного дифференцирования функций // Наука. Инновации. Технологии: научный журнал Северо-Кавказского федерального университета. 2013. №3. С. 7-18.

6. Ланцош К. Практические методы прикладного анализа. М.: Физматлит, 1961. 524 с.