Практикум по матлабу. практикум по матлабу. Физических процессов с использованием
 Скачать 1.13 Mb.
|
|
’PropertyName’ и какие значения оно (свойство) может принимать. Например, в приведенном далее примере, цвет (зеленый -’g’) отрисо- вываемого объекта line c дескриптором h задается оператором set(h, ’color’, ’g’). Для получения значения какого-либо свойства или его присвоения другой перемен- ной используется уже знакомый нам оператор get, использовавшийся в предыду- щем параграфе. Полезно будет заметить, что свойства родительских объектов наследуются объек- тами-детьми. Поэтому при неизбежно поверхностном знакомстве с системой, вно- ся изменения в свойства каких-либо объектов, мы можем получить иногда непред- виденные изменения в работе программ. (Фактически это свойство унаследовано системой MATLAB от языка C++, на котором она написана.) 2.4.2. Отрисовка движущихся кривых MATLAB предоставляет различные способы создания движущихся графиков или анимации. Использование свойства графических объектов ’EraseMode’ (ре- жим стирания) удобно для длинных последовательностей простых графиков, у ко- торых изменение от кадра к кадру минимально. Рассмотрим это на примере дина- мического вывода волны (программа WAVEPAK). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Программа WAVEPAK предназначенная для вывода % % волнового пакета с анимацией изображения % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t=0:0.1:100; % Задание вектора времени x=0:0.3:30; % Задание вектора координат k=1.3; w=0.9; n=length(t); % Вычисление формы волны для момента времени t(1) % во всех точках x y=cos(k*x-w*t(1))+cos(x-t(1)); h=line(x,y); % Подготовка графика волны и присвоение % дескриптора этой линии переменной h % Задание цвета объекта (линии) с дескриптором h % через свойство ’color’ set(h,’color’,’g’); axis([0 30 -3 3]); % Выбор осей координат и их масштаб axis manual; 20 % Сохранение дескриптора линии волны и выбор значения % свойства ’EraseMode’ равным ’xor’. Это сообщает % графической системе MATLAB, чтобы она не % перерисовывала весь график (включая оси, цвет фона % и те точки, которые не изменились), а перерисовала % только те точки, которые изменили свои координаты. set(h,’EraseMode’,’xor’); pause; % Пауза перед запуском движения волны. % Она гарантирует обязательный вывод на экран % накопленного в буфере графики. Дальнейшее % продвижение осуществляется нажатием на любую клавишу. % Можно задать pause(0), тогда нажимать клавишу не надо % Основной цикл вычисления и вывода движущейся волны for i=2:n; % Вычисление формы волны для момента времени t(i) y=cos(k*x-w*t(i))+cos(x-t(i)); % Обновление координат линии, изображающей амплитуду % волны с помощью обновления у объекта (линии) % с дескриптором hсвойства ’XData’ (замена % x-координат на новые значения) и аналогичные % действия с y-координатами. Это стандартный % способ обновления координат точек для анимации. set(h,’XData’, x,’YData’,y); end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Аналогичным образом можно создавать любые движущиеся графики. Пробле- ма возникает тогда, когда скорость компьютера такова, что лимитирующим факто- ром плавного вывода на экран является скорость расчета новых точек. В таких слу- чаях либо улучшают (ускоряют) алгоритм расчета, либо готовят анимацию путем последовательной покадровой записи результатов вывода на экран с последующим прокручиванием получаемого файла в режиме анимации (мультфильма). Задание 5. Попробуйте, с помощью описанного выше алгоритма соединить ани- мацию движущейся волны и графический интерфейс пользователя, описан- ный в п. 2.3 . При этом следует иметь в виду, что необходимо разбить опи- санную выше программу WAVEPAK на две части, первая из которых войдет в основную программу, вызывающую все остальные (типа Drive_Beats), а 21 Рис. 5. Силы, действующие на математический маятник в поле тяжести вторая часть войдет в функцию типа Run_Beats. Если у Вас не получится выполнить это задание самостоятельно, ознакомьтесь с текстами функций в директории WAVEPAK. 3. Об использовании безразмерных переменных Во многих задачах бывает полезно вводить безразмерные переменные. (Это отно- сится, разумеется, не только к работам с применением компьютера.) Показать, как это можно сделать и чем это удобно, лучше всего на конкретном примере. Рассмотрим задачу о движении математического маятника – грузика массы m, подвешенного в поле тяжести на невесомом стержне длины l (рис. 5 ). Будем счи- тать также, что маятник движется в одной плоскости. Пусть на маятник действуют сила трения, пропорциональная скорости грузи- ка, Fтр = −Av, и внешняя переменная сила F (t) = F cos Ωt, направленная горизонтально. (Можно представлять, что маятник помещен в большой плоский конденсатор с вертикальными пластинами, к которым приложено переменное на- пряжение, а на грузике находится постоянный заряд). Для удобства сопоставления с текстом программы мы обозначаем угол отклонения маятника от вертикального направления x. Для угла отклонения нити от вертикали можно записать уравнение ml d 2 x dt 2 = −mg sin x − Al dx dt + F cos x cos Ωt. (1) 22 В отсутствие силы трения и внешней переменной силы период малых колебаний маятника равен, как известно, 2π l/g. Введем вместо времени t новую переменную τ согласно соотношению t = l/gτ ; переменная τ оказывается, очевидно, безразмерной. Тогда уравнение ( 1 ) приводит- ся к виду d 2 x dτ 2 = − sin x − a dx dτ + f cos x cos ωτ, (2) где a = A m l g , f = F mg , ω = Ω l g – (3) безразмерные величины. Такое преобразование обнаруживает некоторые законы подобия: зависимость x (τ) (при заданных x(0), dx/dτ(0)) оказывается одной и той же при разных зна- чениях m, l, g, A, F, Ω, если одинаковы составленные из них безразмерные ком- бинации a, f, ω. Этот факт позволяет существенно сократить объем полного ис- следования задачи, так как достаточно рассматривать различные значения трех па- раметров вместо шести. Иначе говоря, результаты исследования одного маятника можно перенести на другие простым изменением масштабов. Кроме того, при численном определении решения уравнения ( 2 ) мы не будем, как правило, иметь дело с величинами, отличающимися друг от друга на много порядков, в то время как для уравнения ( 1 ) это вполне могло бы получиться при неудачном выборе единиц измерения. И хотя наш компьютер работает с «разма- хом» от -308 порядка ( realmin) до 308 (realmax), лучше эту возможность не использовать без необходимости. Это позволит, в частности, не осложнять грубые оценки, выполняемые без компьютера. Переход от уравнения ( 1 ) к ( 2 ) можно оформить и несколько иначе. Вместо обычных единиц измерения длины, массы и времени выберем «естественные» для данного маятника. За единицу длины примем его длину, за единицу массы – его массу, а единицу времени выберем такой, чтобы было g = 1 (т.е. равной l/g). Подставив l = m = g = 1 в ( 1 ), мы получим уравнение вида ( 2 ), в котором величины t, A, F и Ω подразумеваются выраженными в этих «естественных» еди- ницах. Принимая для этих величин обозначения τ, a, f, ω, приходим к ( 2 ). Чтобы получить на этом пути соотношения ( 3 ), нужно построить из l, m, g множители необходимой размерности. Скажем, сила имеет размерность mg, поэтому запишем F = mgf. Это равенство справедливо в «естественных» единицах и справедливо при переходе к любым другим единицам, если считать f безразмерной величиной, поскольку размерности его левой и правой частей тогда одинаковы. Аналогично мо- гут быть получены и остальные соотношения ( 3 ). Именно выбор естественных для 23 задачи единиц измерения обеспечивает отмеченный выше «умеренный» характер изменения используемых безразмерных переменных. Описанный подход применя- ется весьма часто. В задаче о движении частицы массы m в поле U = −α/R речь может идти как о движении планеты вокруг Солнца, так и о движении электрона вокруг атомного ядра. (В первом случае α = γ m M, где γ – гравитационная постоянная, m – масса планеты, M – масса Солнца; во втором m – масса электрона, α = |qQ|, где q – заряд электрона, |Q| – заряд ядра 8 .) Уравнение движения m d 2 R dt 2 = − α R R 3 . Введем в качестве единицы длины характерную длину R 0 (для астрономической задачи, например, R 0 = 10 8 км, для атомной R 0 = 10 −8 см). Тогда естественно в качестве единицы времени выбрать t 0 = R 3 0 m/α (для движения планеты t 0 окажется порядка года, для движения электрона – порядка периода обращения электрона в атоме). Безразмерные длина r и время τ определяются равенствами t = t 0 τ , R = R 0 r, а уравнение движения приобретает вид 9 d 2 r dτ 2 = − r r 3 . Еще один пример – задачи релятивистской физики частиц, где обычно пола- гают скорость света c = 1. При этом скорость частицы становится безразмерной величиной, равной v/c. Масса же остается размерной величиной, но ее размер- ность не отличается от размерности энергии. Например, масса электрона m = 511 кэВ. 4. Маятник Задачи о колебаниях встречаются во всех областях физики. Во многом колебания совершенно различных физических объектов сходны друг с другом. Простейшие примеры – малые колебания маятника и электрические колебания в цепи, состав- ленной из конденсатора и катушки. 8 Хотя движение электронов в атоме подчиняется не законам движения планет, а законам квантовой механики, необходимо знать также, что дали бы для такого движения законы клас- сической механики. 9 В работе «ПЛАНЕТА» сохранен коэффициент α=1, поскольку предполагается модифи- кация поля U . 24 Такое движение маятника хорошо известно – это гармонические колебания. Закон движения можно записать в виде x = a cos(ω 0 t + ϕ 0 ), где ω 0 – часто- та колебаний, a – амплитуда, ϕ 0 – начальная фаза. (Угол отклонения маятника мы обозначили здесь x. Далее будем использовать безразмерные переменные, о которых говорилось в гл. 3 , т.е. примем массу грузика m и длину маятника l рав- ными единице, а также будем считать ω 0 = 1.) Малые колебания описываются уравнением ¨x + x = 0, линейным относительно функции x, поэтому их обычно называют линейными. В этой работе мы будем исследовать движение математического маятника при больших углах отклонения, специально обращая внимание на отличие от законов малых колебаний. Подобные отклонения называют нелинейными эффектами. 4.1. Свободные колебания Результаты исследования движения маятника удобно представить в виде набора кривых на плоскости (x, p), где p = ˙x – скорость изменения угла. Плоскость (x, p) называется фазовой плоскостью, переменная p – импульсом 10 , а кривые, определяемые параметрически законом движения как x = x(t), p = p(t), – фазовыми траекториями. Фазовая траектория определяется, например, начальными значениями коорди- наты x (0) и импульса p(0). Фазовые траектории линейного осциллятора представляют собой эллипсы, за- даваемые законом сохранения энергии. Для математического маятника это спра- ведливо при малых углах отклонения, в общем же случае, при больших значениях углов отклонения движение математического маятника будет более сложным. Кро- ме колебаний возможно вращение маятника в ту или другую сторону. Угол отклонения маятника достаточно задавать в некоторых конечных пределах, например, принимая −π ≤ x < π. При этом следует представлять, что точки фазовой плоскости (−π, p) и (π, p) отождествлены, иначе говоря, прямые x = π и x = −π склеены друг с другом так, что из полосы −π ≤ x < π, −∞ < p < ∞ 10 Точнее говоря, каноническим импульсом в случае математического маятника называется величина ml 2 ˙x, однако мы учитываем условия m = l = 1. 25 получился цилиндр. Если маятник делает один или несколько оборотов, то точка, изображающая его состояние, движется по кривой, обвивающей этот цилиндр. Отличие колебаний с большой амплитудой от малых колебаний сводится к тому, что закон изменения угла со временем отличен от гармонического, а частота их зависит от амплитуды. При некоторой энергии колебания сменяются вращением. Фазовая траектория, разделяющая на фазовой плоскости области, отвечающие колебаниям и вращению, называется сепаратриссой. (Точнее, это совокупность трех фазовых траекторий, из которых одна – просто точка.) Мы будем сразу же использовать уравнения движения маятника, приведенные к безразмерному виду ( 3 1 ),( 3 2 ). Свободное движение математического маятника без трения описывается дифференциальным уравнением ¨x + sin x = 0. (1) Имеется аналитическое решение этого уравнения (довольно сложное), однако мы будем исследовать движение маятника численно. Запишем это дифференциальное уравнение в виде системы уравнений первого порядка ˙x = p, ˙p = − sin x. (2) Основная идея численного расчета чрезвычайно проста: зная значения коорди- наты и скорости в момент времени t, можно приближенно найти их значения через малый промежуток ∆t x (t + ∆t) = x(t) + ˙x(t)∆t p (t + ∆t) = p(t) + ˙p(t)∆t, взяв значения величин ˙x и ˙p из уравнений ( 2 ). Многократно повторяя такие вычис- ления, мы найд`ем зависимости x (t) и p(t). Б`ольшей точности можно достигнуть, практически не усложняя расчеты, если использовать так называемую вычисли- тельную схему с перешагиванием. В этой схеме вычисляются значения координат в моменты времени t − 1 2 ∆t, t + 1 2 ∆t, t + 3 2 ∆t, t + 5 2 ∆t, ..., а значения скоростей в моменты времени t, t + ∆t, t + 2∆t, t + 3∆t, ... 26 (подробнее об этом сказано в Приложении B ). Данная вычислительная схема и использована в предлагаемой ниже программе. Чтобы отклонение от точного решения уравнения ( 1 ) было небольшим, должен быть достаточно малым шаг по времени ∆t. В любом случае необходим контроль правильности счета (т.е. того, что вычисленные зависимости x (t) и p(t) не слиш- ком сильно отличаются от точных). Одним из методов контроля счета является контроль постоянства полной энергии системы (разумеется, если есть основания считать, что энергия сохраняется). Могут быть и другие сохраняющиеся величи- ны, но часто вообще нет никаких интегралов движения (например, при наличии силы трения). В подобных случаях надежность счета можно проверить, применяя метод повторного счета с уменьшенным шагом ∆t. Если при этом решение остает- ся прежним, значит, шаг был выбран достаточно малым и можно быть уверенным в правильности результата. Далее приводится текст простейшей программы на языке MATLAB, которая является исходной для начала работы над задачей МАЯТНИК. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Фазовая траектория математического маятника % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; % Очистка рабочей области % Задание начальных значений x1=2.2; % Начальная координата p1=0.0; % Начальный импульс dt=0.025; % Шаг по времени axis([-pi pi -pi pi]); % Задание диапазона осей hl=line(x1,p1); % Задание дескриптора линии % Задание параметров выводимой линии % ’EraseMode’,’none’ - это режим вывода без стирания % ’LineStyle’,’:’ - вывод линии в виде пунктира % по умолчанию выводится сплошная линия % ’Color’,’r’ - задание цвета линии set(hl,’EraseMode’,’none’,’LineStyle’,’:’,’Color’,’r’); grid on; % Задание вывода координатной сетки pause; % Пауза, обеспечивает немедленный вывод % рисунка на экран для продолжения программы % необходимо нажать любую клавишу while 1 % Бесконечный цикл 27 |