Brian_-_Matlab_R2007_s_nulya_33 (771739), страница 39
Текст из файла (страница 39)
мдтыв достаточно больп1ое ди, то выбранное значение следовало бы уменышпь, но на самом деле это только ухудшило бы положение. В конечном счете, единг-гвенным параметром в итерации, который мы исполыуехь является постоянная в, а единст асиным недостатком выполнения всех вычислений в М-файле„как мы делали вьппе, является то, что мы не можем видеть промежуточные значения, которые получал~тел в ходе вычисления.
В этом случае мы можем лсуко подсчитать, что в 225 Глава 10. Прикладные задачи подойдет любое значение в между 0 и О.б, и предлагает, что лучшее значение ве, которое следует использовать для данного бх — это значение, которое делает в 0.25. Обратите внимание, что в то время как мы можем получить более точные результаты при уменьшении ))х, если мы уменьшим его в 10 раз, то для компенсации мы должны уменьшить ))К в 100 раз. В этом случае вычисления потребуют в 1000 раз больше времени и в 1000 раз больше памяти! Случай переменной проводимости Ранее мы упоминали, что задача, которую мы решили численно, могла также быть решена аналитически.
Ценность численного метода заключается в том, что его можно применить к подобным дифференциальным уравнениям в частных производных, для которых точное решение найти невозможно или, по крайней мере, оно не известно. В качестве примера мы рассмотрим одномерное тепловое уравнение с переменным коэффициентом, которое представляет неоднородный материал с переменной удельной теплопроводностью К (х) . ди д ( ди) д'и , ди — = — ~~(.) — ~= (.) —,.~(.)— д( дх ~ дх,) дх' дх Для первых производных с правой стороны мы используем симметрическую конечно-разностную аппроксимацию так, чтобы наше дискретное приближение к дифференциальным уравнениям в частных производных станет аМ а л ~ и к и, — и, иу„— 2иу+и,, А;.и — Ус), и „вЂ” и,, Л( Лх' 2Ат 2Ьх гдеК) = К (в + 5))х).Тогда шаг временидляэтогометодабудет и,"." = з)( и,'.'„+ и,"., )+(1 — 2лу )и,".
+ 0.25(л„, — л,, )(и,".„— и,"., ), где в~ = Каза/())х)'. В следующем М-файле, который мы назвали ива(чезп, мы изменим наш предыдущий М-файл, чтобы включить эту итерацию. » курв Кеввчс йцпссйоп ц = Ьеасус(К, с, х, 1псс, Ьс)гу) Ъ Решает одномерное тепловое уравнение с переменным % коэффициентом К на прямоугольнике, описанном векторами х и з с ц(С(1), х] = 1п1С и граничными условиями Дирихле Ъ ц(с, х(1)) = Ьс)гу(1), п(с., х(епс))) = Ьс)гу(2).
1епдСЬ(х); )) = 1епдсЬ (с); Глава 10. Прикладные задачи 227 В этом случае предельное температурное распределение нелинейно. Оно имеет более крутой температурный градиент в середине, где удельная теплопроводность ниже. Напомним, что можно найти точную форму этого предельного распределения, и(х, к) = 20 (1 + (1/и) атссви (х/5) ), приравняв производную к к нулю в исходном уравнении и решив получившееся обыкновенное дифференциальное уравнение. Вы также можете использовать метод конечных разностей, чтобы решить тепловое уравнение в двух или трех координатах. Для этого и других дифференциальных уравнений в частных производных, в которые входят время и две плоскости координат, вы можете также использовать набор инструментов РОВ (ДУЧП), который предоставляет более совершенный «метод конечных элементов». Решение с помошью программы Ь)гпи!(и)( Мы можем также решить тепловое уравнение, используя программу 8)пш!!п(с Чтобы сделать это, мы продолжим приближать производные х с помощью конечных приращений, но будем представлять уравнение как векторнозначное обыкновенное дифференциальное уравнение, где к выступает как независимая переменная.
Программа 8!пш(!пЕ решает модель, используя средство решения ОДУ среды МАТ1АВ, — обе45. Чтобы показать, как это можно сделать, давайте рассмотрим пример, с которого мы начали. В этом случае )с 2 на промежутке -5 5 х 5 5 от времени 0 до времени 4, используются граничные температуры 15 и 25, а начальное распределение температуры 15 — для х < 0 и 25 для х > О. Мы заменяем функцию и(х, с) для постоянного К вектором и, содержащим значения и(х, с), например, для х -5: 5.
Здесь получается 11 значений х, для которых мы отберем и, но так как функция и(х, к) предопределена в конечных точках, мы можем сделать и девятимерным вектором, и в конце мы только присоединим значения в конечных точках. Так как мы заменяем отношение азгйх' его конечно-разностной аппроксимацией и для простоты мы взяли Ох = 1, то наше уравнение становится векторнозначным ОДУ ди — = /с(Аи+с) д( Здесь правая сторона представляет наше приближение к 1(д и I дх ). Матрица 2 2 а выглядит следующим образом -21...0 1 -2 МАТСОВ 228 так как мы заменяем д иl дх в (п, к) наи(п - 1, к) - 2и(п, с) + и(п + 2 2 1, С).
Мы представляем эту матрицу в обозначениях среды МАТ(АВ как 2*еуе (9) + Жар( (81), 1) + 41ад( (8, 1), -1). Вектор о получается из граничных условий, и его первый элемент равен 15, а последний элемент равен 25, элементы между ними — О. Мы представим это в обозначениях среды МАТ(АВ как (15) иегов (7,1) ) 25]. Формула для с получается, исходя из того, что и(1) 2 2 представляети(-4, к) и д и/дх вэтойточке приближается к и(-5, С) — 2и(-4, К) + и(-3, К) = 15 — 2и(1) + и(2) также и в другой конечной точке. Ниже показана модель программы З)пщОп)с, представляющая это уравнение: » орех вувке3в Ьеаке((.вк]1 Обратите внимание, что необходимо определить начальные условия для и как свойства блока (п(вцга!ог (Интегратор), и что в диалоге В)ос)г Рвгвте(вгв (Свойства блока) для блока Ов)п (Приращение) требуется установить вид умножения Мв(пх (Матричное).
Значения с и (1) по и (4) представляют решение и (х, К) в х -4 до -1, азначенияси(б) дои (9) представляют решенияи (х, С) вх = 1до 4, мы принимаем начальное значение и равным (15эопев (4, 1) к 20) 25*опев (4, 1) ]. (Мы используем число 20 как компромисс в х = О, так как это значение справедливо в середине областей, где и равно 15 и 25.) Результат выполнения модели отображается в блоке Зсорв (Индикатор) в виде графиков различных элементов и как функции от с, но более полезно сохранить результат в рабочем пространстве среды МАТ( АВ и затем построить его с помощью команды вцг(.
Кстати, это помогает переустанавливать пределы оси 2 на блоке Зсоре (Индикатор), чтобы переходить от 15 к 25. Чтобы внести эти изменения и запустить модель от в = 0 к С = 4, мы выполним следующие команды » вес рагавг( 'Ьеавец/Бсоре', ' зк(1п', '15' ) ) » еек рагавг('Ьеавец/Зсоре', '"ХВ(ах', '25')к » (соис, иоис] в1гв('Ьеасе<~.вк]1', 10 4])) » в1хр1ов(коис, иоик) МАт).АВ гзо Решение с помошью команды руре Среда МАТ(.АВ содержит встроенную команду рйере для решения дифференциальных уравнений в частных производных в одной координатной плоскости (содержащие и время Ф).
Чтобы узнать больше об этой функции, прочитайте оперативную справку, посвященную рйере. Инструкции по использованию рнера весьма подробны, но требуют математической подготовки. Метод, который использует данная команда, в какой-то степени похож на метод, использованный нами ранее в решении с помощью программы 8(ти)(пК То есть команда использует средство решения ОДУ в е и в конечных приращениях х. Следующий М-файл решает вторую задачу, рассмотренную выше, с переменной проводимостью. Обратите внимание па использование описания функции и подфупкций. » суре Ьеавецехг йипссьоп Ьеасес(ехг Ъ Решает задачу Дирихле для теплового уравнения в стержне, % теперь с переменной проводимостью, 21 узел сетки.
а1 = 0; Это просто говорит, что геометрия линейна. х = 11пзрасе(-5,5,21); = 11пзрасе(0,4,81); во1 = рйере(ш,арйех,$рйехьс,9рйехЬс,х,с); Ъ Выделяет первый элемент решения как и. и = во1(:,:,1); Ъ Объемное отображение часто помогает найти решение.
вигй(х,п,и) сгс1е('йишегьса1 зо1исьоп сошрисей иьсЬ 21 шевЬ роьпсв 1п х. ') х1аЬе1('х'), у1аЬе1('Г'), к1аЬе1('и') йипссьоп с=1; й= (1+ в=О; Ъ [с, й, в] = рйех [х, С, и, Ои()х) (х/5)."2)*Эи()х( йипссгоп иО = рйехгс(х) Ъ Ъ начальное условие при с = 0 иО = 20+5*вйдп(х); йипспйоп (р1,с(1,рг,с(г) = рйехЬс(х1,и1,хг,иг,С) Ъ значения с) равны нулю из-за условий Дирихле е р1 = 0 — в левой конечной точке, рг = 0 — в правой точке р1 = и1-15; Обратите внимание, насколько этот график похож на изображение, полученное на- ми ранее для постоянной проводимости )с ю 2. Мы оставляем читателю возмож- ность поэкспериментировать с моделью для случая переменной проводимости.
Глава 10. Прикладные задачи рл = их-25 с1г = 0; После выполнения мы получим » Ьеаеецех2 с',"„о;" "' "'"„;,„";„; " "-' ""';,;;"„'' -',," о ' '"сз'"-' 'с;. е ' ""' ' о::. ~ . '-"' " "' ' '7 „..".;,'л„, и,::.= "б" ':: .:, ':,".'::'л;:„':.; ~ О. мдт.дв 232 (*) где Š— это время,2 — время реакции,ц — это положение л-ого автомобиля, а «коэффициент чувствительности» Л может зависеть от ц г(е) - ц„(с) (от промежутка между автомобилями) и/или от 6„(е) (от скорости и-го автомобиля). Смысл этого уравнения заключается в следующем.