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

Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. — 2009. — № 1 (18). — С. 122-132

Математическое моделирование

УДК 536.2(075)

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

Е. В. Стефанюк, И. В. Кудинов, Е. В. Ларгина

Самарский государственный технический университет, 443100, Самара, ул. Молодогвардейская, 244.

E-mail: stef-kateSyandex.ru

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

Ключевые слова: нелинейные уравнения, дополнительные граничные условия, приближённые аналитические решения.

Получение решений нелинейных обыкновенных дифференциальных уравнений первого порядка в случаях, когда они включают члены, содержащие искомую функцию в степенях n ^ 3, представляет значительные трудности не только для аналитических, но и численных методов. При достаточно больших n, когда нахождение аналитических решений вообще не представляется возможным, весьма затруднительным оказывается и получение их численных решений. Это связано с тем, что в результате конечно-разностной аппроксимации исходного дифференциального уравнения относительно искомой функции получаются большие системы алгебраических полиномов высокого порядка, методы решения которых (как численные, так и аналитические) пока ещё не разработаны.

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

Стефанюк Екатерина Васильевна — ассистент кафедры теоретических основ теплотехники и гидромеханики; к.т.н. Кудинов Игорь Васильевич — студент.

Ларгина Евгения Валериевна — ассистент кафедры теоретических основ теплотехники и гидромеханики.

порядка. С целью оценки точности предлагаемой методики найдём сначала решения обыкновенных дифференциальных уравнений, имеющих точные аналитические решения. И, в частности, рассмотрим краевую задачу о движении фронта температурного возмущения, возникающую как промежуточный этап решения нестационарной задачи теплопроводности для бесконечной пластины при симметричных граничных условиях первого рода [1, 2]. Задача о движении фронта температурного возмущения в данном случае содержит нелинейное обыкновенное дифференциальное уравнение первого порядка и начальное условие, характеризующее положение фронта в начальный момент времени:

21 dq(Fo) 245

128 dFo 64q(Fo)

(0 < Fo < Foi;0 < q(Fo) < 1); q(0) = 0,

(1) (2)

где q(Fo) —функция, характеризующая перемещение фронта температурного возмущения по безразмерной пространственной координате £ = ^ (0 ^ £ ^ 1) во времени Fo = jp (рис. 1); Fo — число Фурье; a — коэффициент температуропроводности; т — время; R — половина толщины пластины; x — координата; Foi — время окончания первой стадии процесса — стадии нерегулярного режима теплообмена. На рис. 1: 9 = т~-Т0—относительная избыточная температура, где То — начальная температура, Тст — тем-

1,0 е

\ \F° < то) > = Foi

£

R

1,0

Рис. 1. Схема перемещения фронта температурного возмущения по координате £ во времени Ро.

пература стенки при £ = 0, Т — температура.

Уравнение (1) получено при нахождении решения нестационарной

задачи теплопроводности в пятом приближении первой стадии процесса. Это решение имеет вид [1, 2]

14

e(£,Fo) = 5>fc (l

k=o

(3)

где Ао = 1; А = -245/64; А2 = 0; А3 = 455/32; А4 = 0; А5 = -3003/64; А6 = 0; А7 = 2145/16; Аз = 0; А9 = -35035/64; А10 = 1001; Аи = -28655/32; А12 = 455; А13 = 8085/64; А14 = 15.

Уравнение (1) не представляет затруднений для непосредственного интегрирования, и точное аналитическое решение задачи (1), (2) имеет вид

q(Fo) =

yïÔ5Fo0'5.

(4)

k

Положив q(Foi) = 1, из соотношения (4) найдём время достижения фронтом температурного возмущения координаты £ = 1: Fo = Foi = 0,0214285714.

Соотношения (3), (4) представляют решение в пятом приближёнии задачи теплопроводности для бесконечной пластины, полученное с помощью интегрального метода теплового баланса при использовании дополнительных граничных условий. Отличие этого решения от точного [4] не превышает 0,004% [1].

Найдём приближённое аналитическое решение задачи (1), (2). С целью упрощения процесса получения решения введем новую независимую переменную г] = . Тогда уравнение (1) примет вид

21 dq(n) 245Fo1 , , ч ч

ШЖ = ёШ <5>

q(0) = 0. (6)

Для уравнения (5) должно быть выполнено также условие

q(1) = 1. (7)

Решение задачи (5), (6) примем в виде

q(n) = bn\ (8)

где b, X — некоторые пока неизвестные постоянные. Отметим, что из условия (7) следует, что b =1.

Для определения коэффициента X потребуем, чтобы соотношение (8) удовлетворяло не уравнению (5), а некоторому осреднённому уравнению, т.е. уравнению (5), умноженному на dn и проинтегрированному в пределах от П = 0 до n = 1:

1

21 dq 245 Foi

jiwsi-lû ТГ = 0- (9)

0

Подставляя (8) в (9), после вычисления интеграла находим

70Foi + 3X - 3 = 0. (10)

Уравнение (10) содержит две неизвестные величины Foi и X. Для их определения используем ещё одно уравнение, получаемое путём введения дополнительного граничного условия. Такое условие можно найти, если потребовать, чтобы соотношение (8) удовлетворяло уравнению (5) в точке n = 1, т. е.

21 128

dq(y) dn

_ 245Foi v=l 64ç(l)

Последнее соотношение с учётом (7) приводится к следующему алгебраическому уравнению:

21 А 245 0. (11)

128 Fo1 64

Из совместного решения уравнений (10) и (11) находим Л = 0,5; Foi = = 0,02142857143. Отметим, что полученное таким путём значение Foi до одиннадцати знаков после запятой совпадает с точным его значением, а значение коэффициента Л полностью совпадает с показателем степени числа Фурье соотношения (4).

Соотношение (8) с учётом найденных значений Foi и Л для независимой переменной Fo принимает вид

q(Fo) = kFo0'5, (12)

где к = -Лтг = 6,831300513.

Foi'

Точное аналитическое решение (4) задачи (1), (2) можно представить в виде формулы (12), где к = |\/Ю5 = 6,831300511. Отсюда следует, что отличие приближённого аналитического решения задачи (1), (2) от точного наблюдается лишь в девятом знаке после запятой в значении коэффициента k.

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

(3qA + 37q3 + 16g2 - 336q)-^- + 420g + 3360 = 0

dFo

(0 < Fo < Foi; 0 < q(Fo) < 1); q(0) = 0.

(13)

(14)

Задача (13), (14) возникает как промежуточный этап решения нестационарной задачи теплопроводности для бесконечного цилиндра при граничных условиях первого рода [1, 2]. Уравнение вида (13) получается при решении указанной задачи во втором приближении первой стадии процесса.

Непосредственное интегрирование уравнения (13) при начальном условии (14) приводит к соотношению

736 ÎÔ5

Ы 1 + | ) -

92q Пд2 105 + TÔ5"

13g 1260

А

560

Fo = 0.

(15)

Полагая д(Ро!) = 1, из (15) находим время окончания первой стадии процесса: Ро1 = 0,042124470264.

Соотношение (15) относительно неизвестной функции д(Ро) является трансцендентным уравнением, аналитическое решение которого не представляется возможным. Численные значения q в зависимости от Ро приведены в табл. 1.

Таблица 1

ч 0,142 0,323 0,462 0,571 0,665 0,749 0,828 0,902 0,972 1,0

Fo 0,001 0,005 0,01 0,015 0,02 0,025 0,03 0,035 0,04 0,04212

Для нахождения приближённого аналитического решения задачи (13), (14) приведём уравнение (13) к безразмерному виду посредством использо-ваиия новой независимой переменной г] = Тогда задача (13), (14) примет

вид

(Зд4 + 37д3 + 16д2 - 336д)$^ + 420^01 + ЗЗбОЕо! = 0 (16)

ап

(0 < п < 1;0 < ^ < 1);

9(0) = 0. (17)

Решение задачи (16), (17) примем в виде

9(п) = ЬпЛ • (18)

Из условия 9(1) = 1 находим Ь =1.

Потребуем, чтобы соотношение (18) удовлетворяло осреднённому уравнению (16). Для этого проинтегрируем уравнение (16) в пределах от п = 0 до П = 1, т. е.

1

dq

dr] = 0. (19)

(3g4 + 37g3 + 16g2 - 336g) ^ + 420gFoi + 3360Foi

dn

0

Подставляя (18) в (19), после вычисления интеграла получаем

9169(1 + X1) - 226800Fo1 - 201600Fo^1 = 0. (20)

Уравнение (20) содержит две неизвестные величины Fo1 и X1. Для их определения необходимо еще одно уравнение. В качестве такого уравнения используем дополнительное граничное условие, согласно которому в точке П = 1 должно выполняться уравнение (16), т.е.

(3g4 + 37q3 + 16g2 - 336g)^ + 420gFoi + 3360Foi

dn

= 0. (21)

n=1

Подставляя (18) в дополнительное граничное условие (21), с учётом того, что q(1) = 1, относительно неизвестных Fo1 и X1 получаем следующее алгебраическое линейное уравнение:

14X1 - 189Fo1 = 0. (22)

Из совместного решения уравнений (20), (22) находим Fo1 = 0,04212447026; X1 = 0,568680348566.

Отметим, что величина Fo1 до двенадцатого знака после запятой совпадает с точным его значением, полученным из соотношения (15). Решение задачи (16), (17) принимает вид

q(n) = nAl. (23)

Переходя к независимой переменной Fo, получаем

q(Fo) = &1FoAl, (24)

где кг = -Д- = 6,056500253089.

Ро11

Максимальное расхождение результатов расчётов по формуле (24) в сравнении с точным решением (15) составляет 2%.

Для повышения точности получаемого решения применим итеративный метод, заключающийся в следующем. Решение уравнения (16) принимается в виде (18) и определяются значения Л2, А3, Л4, ..., Л10 для случаев, когда фронт температурного возмущения находится соответственно в точках 0,9; 0,8; 0,7; ...; 0,1 координаты £. Затем находится средний коэффициент (Аср) из всех полученных. Причём вклад каждого из полученных коэффициентов в определение среднего значения Л принимается различным. Вклад коэффициента А1, найденного для q(Fo) = 1, принимается за единицу, т.к. этот коэффициент найден с учётом всего пути, проходимого фронтом температурного возмущения (0 ^ £ ^ 1). Соответственно вклад коэффициента Л2, найденного для q(Fo) = 0,9 (0 ^ £ ^ 0,9), принимается равным 0,9, коэффициента Аз, найденного для q(Fo) = 0,8 (0 ^ £ ^ 0,8), — равным 0,8 и т.д. В итоге средняя величина коэффициента Л будет определяться по формуле

. _ А1 + 0,9А2 + 0,8Лз + ••• +0,1Лю

Аср — -, (¿о)

ц

где ц = 1 + 0,9 + 0,8 + ■ ■ ■ + 0,1 = 5,5.

Описанный выше итерационный способ получения решения введён авторами без какого-либо строгого математического обоснования. Однако его высокая эффективность показана на достаточно большом числе конкретных примеров решения нелинейных обыкновенных дифференциальных уравнений [1, 2, 5]. Такой результат можно объяснить тем, что в данном случае удается получить решение, которое наилучшим образом удовлетворяет исходному дифференциальному уравнению на всем отрезке изменения независимой переменной 0 ^ £ ^ 1, определяя при этом всего два неизвестных коэффициента Аср и к. Первый из этих коэффициентов характеризует кривизну кривой q(Fo), а второй — её наклон.

Для определения коэффициента Л2 применительно к уравнению (13) вводится новая независимая переменная г] = щ-, где Ес>2—время, за которое фронт температурного возмущения достигает координаты £ = 0,9. Математическая постановка задачи в данном случае имеет вид

(3<?4 + 37<?3 + 16<?2 - 336<?)$^ + 420<?Ео2 + 3360Ро2 = 0 (26)

ап

(0 ^ п < 1;0 < q(n) < 0,9);

q(0) = 0. (27)

Решение задачи (26), (27) принимается в виде

q(n) = ЬпЛ. (28)

Из условия q(1) = 0,9 получаем Ь = 0,9.

Потребуем, чтобы соотношение (28) удовлетворяло осредненному уравнению (26):

(3g4 + 37g3 + 16g2 - 336g) ^ + 420gFo2 + 3360Fo2

dn

dn = 0.

(29)

Подставляя (28) в (29), получаем

41,922927(1 + X2) - 2Fo2(623 + 560X2) = 0.

(30)

Уравнение (30) содержит две неизвестные величины Fo2 и Л2. Для нахождения ещё одного уравнения используем дополнительное граничное условие, задаваемое в точке п = 1. И, в частности, потребуем, чтобы в этой точке соотношение (28) удовлетворяло уравнению (26), т.е.

(Зд4 + З7д3 + 16д2 - 336g) ^ + 420gFo2 + 3360Fo2

0.

(31)

П=1

Подставляя (28) в (31), с учётом того, что д(1) = 0,9, относительно неизвестных Fo2 и Л2 получаем алгебраическое линейное уравнение

3738Fo2 - 234,44883X2 = 0.

(32)

Из совместного решения уравнений (30), (32) находим Fo2 = 0,03490821620; Л2 = 0,5565688349.

Аналогичный путь решения для д(п) = 0,8; 0,7; 0,6; ...; 0,1 приводит к следующим значениям коэффициентов Fo и Л:

Fos = 0,02817861687 Fo4 = 0,02201149563 Fo5 = 0,01647827834 Foe = 0,01164582637

X3 = 0,5461082219 X4 = 0,5370446616 X5 = 0,5291777391 Xe = 0,5223475478

Foz = 0,007576300062; X7 = 0,5164253894;

Fos = 0,004327043073; Xg = 0,5113069740; Fog = 0,00195048068; Fo10 = 0,0004940273676;

X9 = 0,5069073780; X10 = 0,5031572611.

Для среднего значения X, определяемого по формуле (25), имеем

Xср = 0,5405000014509.

Приближенное аналитическое решение задачи (13), (14) приводится к виду

q(Fo) = fcFoAcp, (33)

где k = Fo1-1 = 0,04212447026-1 = 5,539099039974.

Результаты расчётов фронта температурного возмущения q(Fo) по формуле (33) даны на рис. 2.

1

f

0,8 0,6 0,4 0,2

1\ ^ g^ÇÎFo)

у?-

Z . яг ...........

/

0 0,01 0,02 0,03 №

Рис. 2. Перемещение фронта температурного возмущения ц во времени Ро: 1 — рассчитанное по формуле (15) (точное решение); 2 — по формуле (24); 3 — по формуле (33)

е

0,6 0,2

-0,2

0 0,01 0,02 0,03 Ро

Рис. 3. Изменение относительной невязки е уравнения (13): 1 — для решения по формуле (24); 2 — по формуле (33)

На рис. 3 представлены результаты расчётов относительной невязки уравнения (13) для случаев использования в качестве приближённых аналитических решений соотношений (24) и (33). Относительная невязка определялась по формуле

ц

где ц = 3360 — свободный член уравнения (13); £ — невязка уравнения (13).

Анализ полученных результатов позволяет сделать вывод, что отличие значений q(Fo), найденных по формулам (24) и (33), от точных их значений, в диапазоне чисел 0,005 ^ Fo ^ Foi не превышает соответственно 2 % и 0,05 %.

Рассмотренный выше метод был применен также к решению следующей краевой задачи, имеющей точное решение:

4 _ 2 2,5

¿l.2i"-m"J = t <34>

q(0) = 0, (35)

где q(x) —перемещение фронта температурного возмущения по безразмер-

1 R

ной поперечной координате У = J? (0 ^ у ^ 1) в зависимости от величины

безразмерной продольной координаты ж = 6рдД (0 ^ х ^ Ж1); £ — поперечная координата; п — продольная координата; Я — половина ширины канала; Ре = —число Пекле; шср—средняя скорость течения жидкости; а — коэффициент температуропроводности жидкости.

Задача (34), (35) возникает как промежуточный этап решения задачи теплообмена при ламинарном течении жидкости в плоском канале [5]. И, в частности, уравнение (34) получено при нахождении решения указанной задачи во втором приближении первой стадии процесса, которое имеет вид

вС-.-) = 1-|| + «(|)4 + |(|)в, ре,

где г = 1-у; 6 =

Интегрируя уравнение (34) при граничном условии (35), получаем

Щ<?4 - |<?3 = -840*. (37)

Положив в (37) q(жl) = 1, находим значение координаты ж1 = 0,0209325396825, при котором фронт температурного возмущения достигает центра канала г = 1 (см. рис. 1, где £ соответствует г, а Fo — переменной ж).

Соотношения (36), (37) определяют решение задачи теплообмена в жидкости во втором приближёнии первой стадии процесса. Неудобство формулы (37) состоит в том, что величина q (также как и в соотношении (15)) не может быть выражена в явном виде. В связи с этим, приходится использовать табличные (численные) значения q от ж (см. табл. 2), что несколько снижает ценность соотношения (36) как аналитического решения.

Таб лица 2

ч 0,02093 0,01 0,005 0,001 0,0005 0,0001

X 1,0 0,7693 0,6039 0,3475 0,2745 0,1594

Используя рассмотренный выше метод, применительно к задаче (34), (35) на первом шаге решения получаем соотношение

q(ж) = к1ЖЛ1, (38)

где к1 = 4,004628198; А1 = 0,3588435374; Ж1 = 0,0209325396825.

Отметим, что величина Ж1 полностью совпадает с точным её значением, полученным по соотношению (37).

Более точный изложенный выше итеративный путь решения приводит к следующей формуле для q(ж):

q(ж) = кжЛср, (39)

где к = 3,872183095; Аср = 0,3501450457.

Анализ полученных данных позволяет сделать вывод, что расхождение результатов расчётов по формуле (38) с точным решением, найденным по формуле (37), составляет менее 1%. Расхождение значений q(ж), найденных по формуле (39), с точным решением находится в пределах 0,1 %.

Для уравнений, не имеющих точных аналитических решений, оценку точности получаемых с помощью рассмотренной выше методики приближённых аналитических решений можно выполнять путём анализа невязки исходного дифференциального уравнения. Расчёты показывают, что относительная невязка уравнения (34) при использовании решения (38) в диапазоне 0,00005 ^ ж ^ ж1 = 0,02093 составляет менее 0,05, а при использовании решения (39) —менее 0,01.

По результатам проведенных исследований можно сделать следующие выводы.

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

2. Удовлетворительная точность получаемых решений (отличие от точных менее 0,1%) объясняется тем, что в итерационном процессе получения решения на каждом шаге итерации благодаря использованию дополнительных граничных условий исходное дифференциальное уравнение выполняется точно во всех принятых точках независимой переменной, число которых равно числу шагов итерационного процесса.

3. Ввиду того, что в рассмотренной методике решение принимается в виде степенной функции q = кИ)Л, где подлежащие определению коэффициенты к и Л характеризуют соответственно наклон кривой изменения искомой функции и её кривизну, данная методика может быть применена лишь для решения нелинейных уравнений с монотонным характером изменения искомой функции. В любых других случаях применение дополнительных граничных условий может быть полезным при условии, что искомое решение будет содержать два или несколько членов ряда [1, 2, 5].

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

1. Кудинов В. А. Аверин Б. В., Стефанюк Е. В. Теплопроводность и термоупругость в многослойных конструкциях: Учебн. пособ. для вузов. — М.: Высшая школа, 2008. — 391 с.

2. Кудинов В. А., Стефанюк В. А. Задачи теплопроводности на основе определения фронта температурного возмущения // Извест. АН. Энергетика., 2008. — №5. — С. 141-157.

3. Гудмен Т. Применение интегральных методов в нелинейных задачах нестационарного теплообмена / В сб.: Проблемы теплообмена. — М.: Атомиздат, 1967. — С. 41-96.

4. Лыков А. В. Теория теплопроводности. — М.: Высшая школа, 1967. — 600 с.

5. Кудинов В. А., Стефанюк Е. В., Антимонов М. С. Аналитические решения задач теплообмена при течении жидкостей в плоскопараллельных каналах на основе определения фронта температурного возмущения // Инженерно-физический журнал, 2007. — Т. 80, №4. — С. 176-186.

Поступила в редакцию 28/1/2009; в окончательном варианте — 03/111/2009.

Stefanjuk E. V., Kudinov I. V., Largina E. V.

MSC: 80A17, 80M25

CONSTRUCTION OF THE APPROXIMATED ANALYTICAL SOLUTIONS FOR NONLINEAR ORDINARY DIFFERENTIAL EQUATIONS BASED ON APPLICATION OF ADDITIONAL BOUNDARY CONDITIONS

E. V. Stefanjuk, I. V. Kudinov, E. V. Largina

Samara State Technical University,

244, Molodogvardeyskaya str., Samara, 443100.

E-mail: stef-kateSyandex.ru

Procedure of obtaining of the approximated analytical solutions for nonlinear ordinary differential equations of the first order is developed based on application of additional boundary conditions. The physical value of application of additional boundary conditions is found in solving the initial differential equation in the boundary points of a problem. Application of iterative method allows to solve this equation in the whole range of independent variable meanings and to receive analytical solutions of satisfactory accuracy by averaging the required unknown factors.

Key words: nonlinear equations, additional boundary conditions, approximated analytical solutions.

Original article submitted 28/I/2009; revision submitted 03/III/2009.

Stefanjuk Ekaterina Vasil'evna, Ph. D. (Tech.), Dept. of the Theoretical Basics of Heattech-nics and Hydromechanics.

Kudinov Igor' Vasilievich, Student, Dept. of the Theoretical Basics of Heattechnics and Hydromechanics.

Largina Eugeniya Valerievna, Assist., Dept. of the Theoretical Basics of Heattechnics and Hydromechanics.