Brian_-_Matlab_R2007_s_nulya_33 (771739), страница 38
Текст из файла (страница 38)
Маятник колеблется назад и вперед, но с каждым колебанием угол отклонения становится меньше, пока маятник совсем не вернется в состояние покоя ко времени 20. Одновременно скорость также периодически изменяется, достигая своих наибольших значений в течение каждого колебания, когда маятник находится в середине своего колебания (угол близок к нулю), и достигает нуля, когда маятник находится в конце своего колебания. Увеличение начальной скорости Мы увеличим начальную скорость до 10 . » [в, ха] = ойе45(сс, [Ос0.01с201, [О 101)) » р1ок(ха(с,1), зса(с,2)) матов 218 ш 3 ш 16 На этот раз угол увеличивается до значения, превышающего 14 радиан, перед тем, как кривая подходит к точке примерно (12.5; О). Если быть более точными, кривая закручивается к точке (4п; О), потому что 4п радиан представляет то же самое положение маятника, что и 0 радиан.
Маятник выполнил два полных оборота перед началом своих затухающих колебаний. Скорость сначала уменьшается, но затем повышается после того, как угол проходит через и, поскольку маят:- ник проходит вертикальное положение и получает импульс. Импульса маятника будет достаточно только чтобы еще раз пройти через вертикальное положение в угле Зп. Нахождение начальной скорости, которая заставляет маятник совершать полные врашения Предположим, что мы хотим с точностью до 0.1 найти наименьшее значение начальной скорости, которая требуется, чтобы заставить маятник, начинающий движение из своего исходного положения, выполнить полное вращение один раз.
Будет полезно показать решения, которые соответствуют нескольким различным начальным скоростям на одном графике. Сначала мы рассмотрим целые значения скорости в промежутке от б до 10. » )зо1а » кок а ю Вг10 [В, ха) ю оаа45(д, [0~0.01:20), [О а)) г р1оС(ха(з, 1), зса(1, 2)) » )зо1а оЫ Глава 10. Прикладные задачи 219 1О -2 -6 О 6 1О 16 Начальные скорости 5, б и 7 не являются достаточно большими, чтобы увеличить угол более и, но начальные скорости 8, 9 и 10 достаточны, чтобы заставить маятник сильно качаться. Давайте посмотрим, что происходит на промежутке между 7 и 8.
Второй шаг » )зоЫ » еок в ~ 7.010.219.0 ха] озта45(9, [010.01120], (О а] ) З р1ок(ха(з, 1), ха(з, 2)) а » )о1а оЫ 1О о 6 Мы видим, что ответ находится где-то между 7.2 и 7А. Давайте выполним еще одноугочнение. Третий шаг » )зо16( » еок а ю 7.260.0567.4 [С, ха] = ос]е45(д, [010.01120], [О а]); р1ое(ха( з, 1), хв( з, 2) ) епс] » )зо1зз окк 220 МАТ).АВ о 5 10 15 Мы пришли к выводу, что наименьшая необходимая скорость находится где-то между 7.25 и 7.3. Ф Численное решение теплового уравнения В этом разделе мы будем использовать среду МАТсАВ для численного решения теплового уравнения (известного также как уравнение диффузии), а также для решения дифференциального уравнения в частных производных, которое описывает много физических явлений, например, проводящий тепловой поток или диффузию примесей в неподвижной жидкости.
Вы можете представить протекание диффузии как каплю краски, распространяющуюся в стакане воды. (До некоторой степени вы могли бы также изобразить сливки в чашке кофе, но в этом случае перемешивание, которое появляется при заливке сливок в кофе, усложняет задачу в силу появления движения жидкости, которое впоследствии может ускориться, если мы станем размешивать кофе.) Краска состоит из большого числа отдельных частиц, каждая из которых неоднократно отскакивает от окружающих водных молекул, следуя в действительности по случайному пути. Так как частиц краски очень много, то их индивидуальные случайные движения образуют чрезвычайно детерминированную общую модель, поскольку краска распространяется равномерно во всех направлениях (здесь мы не учитываем возможный результат от действия силы тяжести).
Подобным образом вы можете представить тепловую энергию, распространяющуюся благодаря случайным взаимодействиям соседних частиц. В трехмерной среде тепловое уравнение выглядит следующим образом Здесь и является функцией от с, ж, у и и. Эта функция представляет температуру, или концентрацию примесей в случае диффузии, во время С и для положения в среде (ж, у, я). Постоянная М зависит от используемых материалов и называется удельной теплопроводностью в случае рассмотрения теплового потока или Глава 10. Прикладные задачи 221 коэффициентом диффузии для задачи, связанной с явлениями диффузии.
Чтобы упростить задачу, давайте предположим, что среда не трехмерна, а одномерна. В этом случае можно рассмотреть диффузию в тонкой заполненной водой трубе или тепловой поток в тонком изолированном проводе. Давайте, прежде всего, рассмотрим случай теплового потока. Тогда дифференциальное уравнение в частных производных примет вид ди ди 2 ' д( дх' гдеп (х, с) -температуравовремя с нарасстоянииж вдольпровода. Решение методом конечных разностей Чтобы решить это дифференциальное уравнение в частных производных, нам необходимы оба начальных условия в форме и (х. 0) = к (ж), где Е (х) дает температурное распределение в проводе во время 0 и граничные условия в конечных точках провода, назовем их х = а их = Ь.
Мы выбрали так называемые граничные условия Дирихле и (а, к) = з., и ц (Ь, с) = кь, которые соответствуют температурам, равным постоянным значениям Т, и Ть в этих двух конечных точках. Хотя в данном случае можно получить и точное решение, позвольте нам показать численный метод конечных разностей. Для начала отметим, что мы можем хранить на компьютере последовательность значений температуры и только в дискретном наборе времен и для дискретного набора положений х. Пусть времена буд)т О, ПК, 2йс ... Ийс, и пусть положения будуте, а + йж ... в + Упж = Ь, ам~ = п (а + зпк, пйс). После переписывания дифференциального уравнения в частных производных в соответствии с условиями конечно-разностной аппроксимации для производных, мы получим 2 Это самые простые аппроксимации, которые мы можем использовать для производных, и этот метод может быть повторно выполнен с использованием более точных приближений, особенно для производной к.
Таким образом, если для отдельного п мы знаем значения и ~ для всех 5, то мы можем решить уравнение, показанное выше, для каждого 5, и,"" = и,". + —,(и,".„— 2и,"+и,"., ) = з(и,".„+ и,"., )+(1 — 2л)и,"., где в = )сбс/(Дж)'. Другими словами, это уравнение говорит нам, как найти температурное распределение в шаге времени и + 1, зная температурное распределение в шаге времени гь (В конечных точках 3 = 0 и з ~ Ю, это уравнение обращается к температурам, находящимся вне ограниченного диапазона для х, но в этих точках мы не будем учитывать данное уравнение, а вместо этого применим МДТ1ДВ 222 граничные условия.) Мы можем воспринимать это уравнение, считая, что темпе.
ратура в данном положении в следующем шаге является средневзвешенным числом его температуры и температур соседних участков в текущем шаге времени. Другими словами, во время йа данный участок провода длины ох передает каждому из своих соседних участков часть своей тепловой энергии в и сохраняет остающуюся часть 1 - 2в. Таким образом, наше численное воплощение теплового уравнения является дискретной версией микроскопического описания диффузии, которую мы привели вначале, когда тепловая энергия распространяется благодаря случайным взаимодействиям между соседними частицами. Следующий М-файл, который мы назвали пеабт, вычисляет описанные выше операции.
» куре Ьеак 1ипсс1оп и = Ьеас(К, с, х, 1п1с, Ьдгу) Ъ Решает одномерное тепловое уравнение на прямоугольнике, Ъ заданным векторами х и С с начальным условием и(С(1),х)= 1п1С $ и граничными условиями Дирихле и (С, х (1)) = Ьс)гу (1), Ъ и (с, х (епс))) = Ьдгу (2). Т = 1епдСЬ(х); ы = 1епцсЬ (с); дх = теап (с)1ЙЙ(х))) бе = шеап (с)1Н (С)); в = К*с)с/с)х"2; и = хегов ()Ч, Т); и (1,: ) = 1п1с бог и = 1:И-1 и(п+1,2:Ю-1) = в*(и(п, 3: Т) + и(п,1:,Т-2) ) + (1-2*в) *и (и, 2: Т-1) ) и (и+1, 1) = Ьс)гу(1); ц (и +1, Т) = Ьс)гу(2) епс) Функция принимает в качестве входных параметров значение К, значения векторов к и ж, вектор начальных значений 1пЫ (предполагается, что его длина равна длине вектора ж), а также вектор Ьбк1г, содержащий пару краевых значений.
В результате выполнения этой функции получается матрица значений и. Обратите внимание на то, что индексы массивов в среде МАТОВ должны начинаться с 1, а не с О. Мы немного отклонились от нашего более раннего примечания, что и 1 Глава 10. Прикпадные задачи 223 соответствует начааьноьпу времени, а 3 = 1 соответствугп левой конечной точке. Обратите янимание так:ке на и, что в первой строке по~шс оператора цикла (ог мы вычисляем весь рял и за исклкшг.нием первых и последних значений, я одной с.грокгс кая'дый элемент является векторг~м ~~иной У - 2, с индексом 3, ) вели «ейгцям на 1 в элементах я (и, 3: Т) и уменьшенным на 1 в элементах и (и, 1 г,Т-2) .