Практикум по матлабу. практикум по матлабу. Физических процессов с использованием
 Скачать 1.13 Mb.
|
|
x1=x2;p1=p2; end; (Здесь шаг по времени ∆t обозначен dt, pi=π.) Строго говоря, следовало бы перед началом основного цикла вставить операцию x=x-p*dt/2, имея в виду, что угол отклонения маятника x и скорость его изменения p в вычислительной схеме с перешагиванием относятся к разным моментам времени. Можно не делать этого, но следует учитывать, что, начальные значения x и p относятся к разным моментам времени. Работа программы понятна из приведенных в ее тексте комментариев. Следует только сделать несколько замечаний. После запуска программа откроет графиче- ское окно и нарисует в нем одно окно с системой координа, осями и координатной сеткой, после чего перейдет в режим ожидания (команда pause). Для продолже- ния расчета необходимо нажать любую клавишу, и программа начнет вывод фа- зовой траектории. Поскольку используется бесконечный цикл, выполнение про- граммы можно прервать клавишей Ctr+C. Для повторного прогона с другими параметрами необходимо ввести соответствующие изменения в текст программы и перезапустить ее 11 Задание 1. Уменьшая шаг dt, убедитесь, что фазовая траектория воспроизво- дится при этом без изменений в течение нескольких периодов (разумеется, с той точностью, какую допускает графическое изображение). Увеличивая шаг, 11 Можно, конечно, уже сейчас, используя опыт предыдущей задачи (Биения) и разрабо- танные m-файлы, создать графический интерфейс и работать далее с ним. Но можно это сде- лать и позже или не делать вовсе. Это зависит от того, что проще – разрабатывать интерфейс или использовать примитивные приемы остановки и перезапуска программы. 28 достигните такой его величины, чтобы появилось явно видимое искажение формы фазовой траектории. Изобразите несколько разных фазовых траекторий. Задавать начальные точки для вывода нескольких фазовых траекторий мож- но либо используя графический интерфейс (см. п. 2.3 ), либо с помощью по- вторного запуска программы, как это описано выше. Постройте фазовые траектории разного типа, отвечающие колебаниям и вра- щению, а также сепаратриссу. Включите в программу дополнительный цикл, который обеспечил бы постро- ение целой серии фазовых траекторий с разными значениями начальной ко- ординаты (и/или импульса) 12 Задание 2. Исследуйте, с какой точностью выполняется при расчете закон со- хранения энергии. Удобно выводить зависимость энергии E от координаты x, чтобы точка (E, x) изображалась под точкой (p, x) (или над ней) в своем окне, а фазовая траектория -в своем. Для изображения двух разных графиков на одном листе можно открыть подрисунки на одном листе (функция subplot, см. п. 2.1 ). При этом для получения двух изменяющихся одновременно гра- фиков можно использовать такой прием. 1. Вставить в начало программы определение двух подокон и определение двух разных дескрипторов для каждой из линий. subplot(2, 1, 1); % определение 1-го подокна axis([-pi pi -pi pi]); % Задание диапазона осей hl=line(x,p); % Задание дескриптора 1-й линии subplot(2, 1, 2); % определение 2-го подокна axis([-pi pi 0.8 1.2]); % Задание диапазона осей he=line(x,E); % Задание дескриптора 2-й линии 2. После соответствующих вычислений новых значений координат, импуль- сов и энергии внутри цикла нанести на каждый из подрисунков свою но- вую точку с помощью, например, такой последовательности операторов. 12 Интересно не писать внешний цикл, а сделать вектора x и p матрицами и написать аналог приведенной выше программы для матриц, выполняя тем самым в одном цикле расчет целого семейства траекторий. 29 ............................. set(hl,’XData’,x,’YData’,p); set(he,’XData’,x,’YData’,E); ............................ Задание 3. Постройте изображение качающегося (в соответствии с расчетом) маятника, а также зависимость x (t). При изображении качающегося маят- ника обратите внимание на различия в задании свойства ’EraseMode’ -’xor’, ’background’ или ’none’. Для построения движущегося отрезка, изображаю- щего маятник, рекомендуется присвоить (с помощью оператора set) свойству ’EraseMode’ значение ’background’ для соответствующего объекта. Изобразите также (в отдельном окне) зависимость x(t), чтобы можно было видеть движение маятника одновременно в трех разных формах. Задание 4. Получите зависимость частоты колебаний маятника от амплитуды. Для этого удобно начальную точку фазовой траектории (x 0 , p 0 ) задавать на оси p = 0 и определять полупериод колебаний T/2, подсчитывая число ша- гов, необходимых для смены знака у импульса p. Для представления зависимости ω (x 0 ) удобно выделить отдельное графиче- ское окно. Задание 5. Включите в программу силу трения, пропорциональную скорости. Коэффициент пропорциональности следует включить в число параметров, до- ступных оперативному изменению (если вы работаете с графическим интер- фейсом). 4.2. Вынужденные колебания Будем рассматривать только колебания маятника под действием гармонической си- лы (и с учетом силы трения), а также будем иметь в виду ту модель, которая опи- сана в ( 1 ). В этом случае безразмерное уравнение движения имеет вид ¨x = − sin x − 2λ · ˙x + f · cos x · sin Ωt. (3) 4.2.1. Переходные колебания Если маятник в начальный момент покоился, а в дальнейшем на него действует периодическая сила, то в течение некоторого времени он раскачивается, происхо- дит нерегулярное движение. При этом возможен не только рост амплитуды, но 30 и ее периодическое нарастание и убывание. Спустя же какое-то время движение становится установившимся, периодическим, причем с частотой внешней силы Ω. Приближение к такому режиму асимптотическое, т.е. происходит, строго говоря, бесконечно долго, тем не менее вполне можно указать интервал времени, спустя который нерегулярные процессы оказываются незаметны (с обусловленной точ- ностью). Обычно для этого нужно примерно такое же время, как для затухания свободных колебаний. Задание 6. Получите различные режимы переходных колебаний – с монотон- ным ростом амплитуды и с ее осцилляциями. При этом предпочтительно вы- брать небольшое значение амплитуды силы f , такое, чтобы даже попав в ре- зонанс, эта сила не приводила маятник во вращение. Коэффициент в силе трения должен быть не слишком велик, чтобы свободные колебания затухали в течение многих периодов. Одна из возможностей наблюдения эффекта – выводить зависимость x (t) в таком масштабе, чтобы график был сильно сжат в направлении оси t 13 4.2.2. Резонанс Если уравнение ( 1 ) для свободных колебаний математического маятника в принци- пе решается аналитически (решение выражается через специальные функции), то уравнение вынужденных колебаний с трением не решается и его можно исследо- вать качественно или численно. До начала численного решения полезно рассмот- реть качественные особенности предполагаемого решения. Рассмотрим случай ма- лых колебаний и установившийся (стационарный) режим. При малых углах от- клонения маятника (x 1) несложно проанализировать зависимость амплитуды установившихся колебаний маятника под действием гармонической силы F (t) = f cos ωt от частоты силы ω 14 Если заменить sin x на x, то получим линейное относительно x уравнение (гар- моническое приближение): ¨x + 2λ ˙x + x = f cos ωt. Его решение удобно искать в виде x = Ae iωt + A ∗ e −iωt , (4) 13 В этом задании предполагается ограничиться исследованием, не затрагивающим всерьез особенности переходных процессов, связанные с нелинейными эффектами. 14 Множитель cos(x) не играет роли в рассматриваемом нами процессе, поэтому мы этой зависимостью пренебрегаем, тем более, что при x 1 cos(x) ≈ 1. 31 представив при этом силу в форме F = (f/2)(e iωt +e −iωt ). Приравняв слагаемые с одинаковой зависимостью от времени, получим A = f 2(1 − ω 2 + 2iλω) , (5) откуда x = a · cos(ωt + δ), где a = 2 · |A| = f/ (1 − ω 2 ) + 4λ 2 ω 2 , δ = arcctg ω 2 − 1 2λω . Учтем теперь следующий (ангармонический) член разложения sin x, по-прежнему считая угол x малым: ¨x − 2λ ˙x + x − x 3 / 6 = f cos ωt. Дополнительные слагаемые, получаемые за счет подстановки ( 4 ) в член x 3 / 6, бу- дут иметь вид x 3 / 6 = (A 3 e 3iωt + 3A 2 A ∗ e iωt + 3AA ∗2 e −iωt + A ∗3 e −3iωt )/6, Слагаемые, пропорциональные e ±3iωt , при- 0.6 0.7 0.8 0.9 1 1.1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 ω a Рис. 6. Резонансная кривая при на- личии гистерезиса ведут к появлению в ( 4 ) малых добавок бо- лее высокого порядка. Это можно предста- вить себе как появление в правой части урав- нения силы, пропорциональной cos 3ωt, а по- скольку их частота далека от резонанса, то и вклад их будет мал. Теперь, как и раньше, собирая и прирав- нивая члены, имеющие одинаковую зависи- мость от времени, например при e iωt , вместо ( 5 ) получаем (1−ω 2 +2iλω +|A| 2 / 2)A = f / 2, откуда получаем уравнение [(1 − ω 2 + |A| 2 / 2) 2 + 4λ 2 ω 2 ]|A| 2 = f 2 / 4. (6) Это уравнение разрешаем относительно ω 2 , после чего зависимость a (ω) легко представить с помощью MATLAB графически (рис 6 ). Аналитическое исследо- вание этой зависимости приведено, например, в [ 7 ]. 32 Зависимость a (ω) представляется существенно разной при различных значе- ниях λ и f . Однако при относительно больших амплитудах сделанное нами при выводе (3) предположение о малости добавок к гармоническому колебанию нару- шается, так что следует ожидать отклонения от этих кривых тех значений ампли- туд, которые мы можем найти в компьютерном эксперименте (и которые являются, разумеется, истинными). Появляются такие участки кривой, которые могут быть достигнуты при непре- рывном изменении частоты либо только путем ее повышения, либо только путем понижения. Для наблюдения такого явления нужно обеспечить, чтобы при изменении часто- ты Ω на малую величину ∆Ω не происходило резких скачков силы f ·cos x·sin Ωt. Тогда переходной процесс будет завершаться быстро и точка x 0 (ω) будет переме- щаться вдоль одной ветви резонансной кривой (разумеется, исключая тот случай, когда скачок на другую ветвь становится неизбежен). Чтобы задать такое плавное изменение силы, нужно избежать скачка фазы си- лы Ωt при изменении частоты Ω → Ω + ∆Ω. Для этого в программе достаточно ввести переменную, равную этой фазе, скажем wt, и на каждом шаге по времени наращивать именно эту фазу: wt=wt+dwt, где dwt=w*dt – приращение фазы. При изменении же частоты dw = ∆Ω следует просто изменять приращение фазы: dwt=dwt+dw*dt 15 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Программа отрисовки зависимости стационарной амплитуды % от частоты вынуждающей силы %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Обратите внимание на то, что lam -- это половина коэффициента % при импульсе в выражении силы трения! clear; f=0.26; % Амплитуда вынуждающей силы lam=0.095; % Половина коэффициента при импульсе % в выражении силы трения! a=0.05:0.001:1.0; aa=a.^2; d=4*lam^4-4*lam^2*(1-0.5.*aa)+f^2./(4*aa); k=find(d>=0); aa=aa(k); 15 Тем самым в качестве фазы будет использована величина t 0 ω (t)dt. 33 a=a(k); d=d(k); w1=sqrt(1-2*lam^2-0.5.*aa-sqrt(d)); w2=sqrt(1-2*lam^2-0.5.*aa+sqrt(d)); h=plot(w1,2*a, w2, 2*a); axis([0.6 1.0 0 2]) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Задание 7. Получите зависимости амплитуды колебаний маятника от частоты внешней силы в случае малых колебаний и в случае, когда наблюдаются "скач- ки"амплитуды. Наложите эту зависимость на "теоретическую", полученную с помощью уравнения ( 6 ). 4.2.3. О случайном движении Упомянем, наконец, об интересном явлении, имеющем место в исследуемой меха- нической системе. При движении вблизи сепаратриссы (в отсутствие трения) ма- ятник много времени проводит близко к вертикальному, неустойчивому положению грузом вверх. В этих условиях даже очень маленькое изменение начальных значе- ний угла или скорости может привести к тому, что он пройдет через вертикальное положение и сделает оборот, вместо того чтобы качнуться обратно. Очень малень- кое отклонение вырастает во много раз! Поскольку расчеты ведутся непременно с округлением, то спустя некоторое время начальные условия "забываются". Движе- ние в ближайшей окрестности сепаратриссы оказывается в определенном смысле случайным. Разумеется, если задать те же начальные значения x и p, при повторном счете зависимости x (t) и p(t) воспроизведутся полностью. Включение же внешней периодической силы может заметно расширить эту об- ласть случайного движения на фазовой плоскости. Однако наблюдать указанное явление непросто, так как непросто отличить ис- тинно случайное движение от регулярного, но кажущегося случайным. Подробнее вопрос о возникновении хаотического движения в детерминированной системе бу- дет затронут в работе “ШАРЫ”. 34 5. Движение частиц в центральном поле Закон движения частицы r(t) в заданном поле U(r) полностью определяется урав- нением движения m ¨r = − ∂U ∂ r = F(r) (1) и начальными условиями r(0) = r 0 , v(0) = v 0 . В некоторых случаях диффе- ренциальные уравнения ( 1 ) можно решить аналитически, в других возможно каче- ственное исследование. Однако в сколько-нибудь сложных полях U (r) приходится применять численные расчеты. 5.1. Траектория финитного движения В данной работе изучается, в первую очередь, движение в центральных полях. Что- бы ориентироваться в ожидаемых результатах, необходимо вспомнить качествен- ные способы исследования движения в центральном поле (см. [ 7 ]). В этом слу- чае сохраняется момент импульса M = m[rv], а траектория лежит в плоскости, перпендикулярной вектору M. Мы ограничиваемся таким случаем, когда частица не удаляется от центра неограниченно (финитное движение) и не падает в центр поля. Тогда, вообще говоря, траектория оказывается расположенной между дву- мя окружностями r min ≤ r ≤ r max и заполняет это кольцо, проходя как угодно близко к любой его точке. Финитное движение в кулоновском поле U 0 = −α/r, в отличие от этого об- щего случая, есть движение по замкнутой траектории – по эллипсу с полуосями: большой a = α/2 | E | и малой b = M/ 2m | E |. В других центральных полях траектории при выбранных наудачу значениях энер- гии и момента импульса окажутся почти замкнутыми лишь после очень большого числа оборотов 16 Самый простой метод численного решения уравнения ( 1 ) (метод Эйлера) состо- ит в том, чтобы, выбрав малый "шаг" ∆t, непосредственно определить приращения координат и скорости равенствами ∆r = v∆t, ∆v = 1 m F(r(t))∆t. (2) 16 Поле U ∝ r 2 , которое в данной работе не рассматривается, обладает такой же особенно- стью, как кулоновское. 35 Многократно повторяя такое вычисление, получим ряд последовательных зна- чений координат и скорости. Гораздо большей точности можно добиться, если вместо ( 2 ) использовать раз- ностную схему с перешагиванием ∆v = 1 m F(r(t + ∆t)), (3) т.е. при вычислении силы использовать новые значения координат в момент време- ни t + ∆t (подробнее об этом см. в Приложении B ). Далее приводится отрывок текста программы, в которой реализован описанный выше алгоритм для вычисления и изображения на экране траектории частицы в кулоновском поле. В программе alpha обозначает α/m, а dt = ∆t, v=(vx,vy), F/m = a =(ax,ay). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Учебная программа для решения задачи Планеты % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; % Очистка рабочей области % Задание начальных значений |