Численные методы решения уравнение маятника

Решение дифференциальных уравнений численными методами

Цель: Используя численные методы решить заданную условием задачу сводящуюся к дифференциальному уравнению либо к их системе. Исследовать фазовое движение маятника в данной задаче (нарисовать фазовую траекторию), построить осциллограмму колебаний. Научится программировать задачи с участием дифференциальных уравнений. Изучить численные методы решения и программирования таких задач.

Анализ условия, решение в общем виде:

Путём подстановки длины нити в уравнение математического маятника получаем следующую систему уравнений:

Выбор численного метода решения и анализ его сходимости:

Данная задача решается известными нам методами: Эйлера и методом Рунге-Кутта. Для начала решим её методом Эйлера. Метод Эйлера характерен своей простотой – это его основное преимущество. Однако у него малая точность (первый порядок точности). В данном методе ошибка «накапливается» потому что приращение значения функции при переходе к следующему шагу заменяется приращением кординаты касательной к кривой решения от предыдущей точки. Метод Эйлера точнее, когда вторая производная функции имеет как можно меньшее значение и «шаги» очень малы.

Алгоритм решения уравнения:

Описание алгоритма: рассмотрев выше составленные уравнения, решаем их по методу Эйлера. Изменяя время на некоторую величину dt, программа выполняет один «шаг» и отправляет все нужные нам значения в графики. Повторяя это, получаем новые значения.

Блок-схема (Для общего случая):

G, l0,a, w,ww, ww0,f, f0,t, dt, f1:extended;

Function func(t, f:extended):extended;

Procedure TForm1.Button1Click(Sender: TObject);

Нелинейный маятник. 1 Безразмерное уравнение движения физического маятника с вязким трением.

Александра Гаевская 5 лет назад Просмотров:

1 Нелинейный маятник. 1 Безразмерное уравнение движения физического маятника с вязким трением. Уравнение движения физического маятника с учётом вязкого трения: I φ + b φ + mga sin(φ) =, (1) где I момент инерции, b коэффициент момента сил вязкого трения, m масса маятника, g- ускорение свободного падения, a расстоянии от точки подвеса до центра масс. Введя частоту ω, перепишем уравнение (1) в виде: φ + b I mga φ + sin(φ) = I ω ω = mga I φ + b I φ + ω sin(φ) = () И приведём полученное уравнение () к безразмерному виду, удобному для анализа: φ + b ω I φ + sin(φ) = ω τ = ω t, d ω dt = d b, β = dτ ω I φ + βφ + sin(φ) = (3) В уравнении (3) производные берутся по безразмерному времени τ = ω t. При этом максимальное значение амплитуды безразмерной скорости φ соответствующей колебательному движению маятника без затухания составляет: φ max =. Отличие решений уравнения (3) от решения нелинейного маятника без затухания определяется только значением безразмерного параметра β, который в дальнейшем предполагается малым: 2 Фазовая диаграмма нелинейного маятника с затуханием. С помощью численного решения уравнения (3) можно определить фазовый портрет нелинейного маятника для заданных начальных условий и параметра затухания β. На рисунке 1 приведён пример фазового портрет для β =.1 и с начальной безразмерной скоростью φ =. На рисунке также показаны красным цветом сепаратрисы нелинейного маятника без затухания ограничивающие область колебательных решений. Рис. 1: Пример фазового портрета для нелинейного маятника с затуханием вычисленный с помощью численного интегрирования уравнения движения (3)

3 3 Закон сохранения энергии для нелинейного маятника с затуханием Интегрирование уравнения (3) по углу отклонения: φ dφ+βφ dφ+sin(φ)dφ = φ dφ +β φ dφ+ sin(φ)dφ = Const φ cos(φ) = β φ dφ + Const (5) Постоянную интегрирования найдём из начальных условий. В начальный момент времени (τ = ) будем считать, что маятник проходит положение равновесия (φ = ), и имеет скорость φ, тогда: Const = φ 1 (6) В момент первого максимального отклонения маятника: cos( ) = β φ dφ + φ 1 β φ dφ = φ + cos() 1 ( β φ dφ = φ ) sin (7) При возвращении к положению равновесия угловая скорость маятника φ 1 из-за трения станет меньше чем начальная φ : 1 1 = φ 1 β φ φ 1 = φ β φ dφ β φ dφ β φ dφ φ dφ В выражения (7) и (8) входят интегралы учитывающие потери энергии на трение. Потери механической энергии связаны с работой сил трения. (8) 4 Связь амплитуды колебаний и угловых скоростей в положении равновесия. Для оценки значения интегралов в уравнении (8) используем следующее приближение которое, обосновывается на рисунке. 3

4 Рис. : Участок траектории маятника на фазовой плоскости Интегралы, входящие в (8) равны площадям заштрихованных на рисунке областей. Примем следующее приблизительное правило: φ dφ A φ (9) φ dφ A φ 1 которое означает, что отношение площади верхней заштрихованной области к площади прямоугольника со сторонами φ и такое же как и отношение нижней заштрихованной области к площади прямоугольника со сторонами φ 1 и. Знак во втором уравнении введён для учёта отрицательного значения скорости φ 1 (заметим, что оба интеграла положительны). Это приблизительное равенство тем лучше будет выполняться, чем меньше параметр безразмерного затухания. Тогда уравнение 4

5 (8) преобразуется: φ 1 А уравнение (7) приводится к виду: ( βa φ = φ sin φ + φ 1φ = φ 4 sin ( или = φ βa(φ φ 1) (1) φ + φ 1 = βa ) ( βa φ = φ 4 sin ) φ 1 φ = 4 sin ( ( sin ) = φ 1φ 4 ) ) (11a) Здесь знак — в правой части учитывает то, что скорости φ и φ 1 всегда должные иметь разные знаки. Другой способ оценки интегралов (9) состоит в том, что угловая скорость φ в подынтегральном выражении выводится через формулы (5) и (6) но с параметром затухания β =. (Метод последовательных приближений.) Как показано в приложении 8.1, это приводит к следующему соотношению между амплитудой и скоростями: sin ( Среднее между (11a) и (11б): ) φ = + φ 1 8 (11б) sin ( ) = (φ φ 1) 16 (11в) Фактически, (11б) это средний квадрат синуса половинного угла амплитуды, (11в) квадрат среднего синуса половинного угла амплитуды, а (11a) это среднее геометрическое. (Имеется в виду усреднение по двум значениям sin( /) получаемым из модели без учёта трения, β =, для скоростей φ и φ 1.) Отметим, что эти выражения позволяют вычислять амплитуду колебания без явного знания значения параметра затухания. Но информация об этом параметре связана с различием угловых скоростей φ и φ 1. На рисунке 3 показаны графики относительной ошибки выражений (11a), (11б) и (11в) для β =.1 в зависимости от амплитуды, полученные на основе численного решения. Видно, что выражение (11в) является 5

6 Рис. 3: Относительная точность формул (11a), (11б) и (11в) наиболее точным для оценки амплитуды колебаний, относительная погрешность не превышает величины.%. Для значений β 7 Если полностью пренебречь затуханием, тогда из (5), (6) и (7): π 1 + dφ = φ ( ( sin ) ) + 9 ( π 1 + ( 4 ) ( 4 ( sin dφ sin ( ) ( sin φ ) = K ( ( sin ) ) ) ) ( ( sin ( ) ) ) ) (13) Через K() обозначена специальная математическая функция — полный эллиптический интеграл первого рода. При переходе от безразмерных величин обратно к размерным получим следующее выражение для периода нелинейных колебаний: T = τ ( sin( 1 τ 1 τ 1 = = T ω π/t π = T )) π ( ) 11 ( ) ) (14) 4 T ( Здесь T период линейных колебаний маятника. В приложении 8.3 показан способ численного вычисления полного эллиптического интеграла K(x ), который используется для теоретического определения значения периода при заданной амплитуде колебаний. 6 Определение периода линейных колебаний. Измерение периода линейных колебаний маятника можно проводить двумя разными способами. В первом случае проводится прямое измерение периода малых колебаний маятника, когда влияние нелинейности пренебрежимо мало. Во втором случае проводится измерение угловой скорости маятника при прохождении через положение равновесия, при движении из начального состояния с углом отклонения равным 18 градусов и нулевой начальной скоростью. Период определяется косвенно через измеренную скорость. Недостаток первого метода состоит в том, что энергия малых колебаний очень мала и поэтому оказывается существенным влияние сил трения разной природы (сухое трение, флуктуации вязкого трения о воздух связанные с потоками воздуха через спицы маятника). Во втором случае энергия маятника значительно больше, хотя потери на 7

8 трения всё равно не равны нулю. Рассмотрим второй метод более подробно. При отсутствии трения закон сохранения энергии связывает угловую скорость ω x и максимальную амплитуду качаний маятника, которая равна 18 градусам: Iωx = mga Обозначения такие же как и в уравнении (1). Выразим отсюда угловую скорость ω x через угловую частоту линейных колебаний маятника ω : mga ω x = I = ω (15) Таким образом, измеряемая угловая скорость получается в два раза больше чем угловая частота линейных колебаний. Период линейных колебаний: T = π = 4π (16) ω ω x Уравнение (16) справедливо при отсутствии трения, в случае вязкого трения можно оценить его вклад в качестве поправки. В приложении 8. выведено выражение (5) в котором учтено влияние трения в рамках первого приближения. Именно это выражение используется в работе при обработке предварительных измерений. 7 Измерения скорости и полупериода колебаний с помощью оптоэлектронного датчика. Оптоэлектронный датчик располагается так, что спица маятника пересекает луч оптопары тогда, когда маятник проходит через положение равновесия маятника. Спица маятника имеет в месте прохождения луча угловую ширину Δφ; моменты пересечения спицей луча происходят так, как показано на рисунке. Здесь отмечены моменты нарушения T 1 i и восстановления T i оптического тракта оптопары (i = 1,, 3..). Зная угловой размер спицы маятника и учитывая, что этот размер достаточно мал (Δφ.47.), можно определять угловую скорость маятника в моменты прохождения положения равновесия: φ i = Δφ, Δt i (17) Δt i = T i T 1 i 8

9 Рис. 4: Зависимость угла поворота маятника от времени синяя кривая. Красные кривые положение граней спицы маятника. Также показаны моменты времени, регистрируемые с помощью оптопары Полупериод колебания: T 1,i = T i + T 1 i T i 1 + T 1 i 1 (18) Эти выражения используются для определения измеряемых в опытах величин: угловой скорости и полупериода колебаний. Период колебаний определяется как сумма двух последовательных полупериодов. Соответствующая этому периоду амплитуда как среднее арифметическое от амплитуд (формула (11в) ) двух последовательных качаний маятника. Если вместо углового размера Δφ в формуле (17) использовать эффективный угловой размер Δφ n : Δφ n = Δφ ω (19) то угловая скорость будет определяться в безразмерных единицах. В опытах с качанием маятника с начальной амплитудой в 18 градусов угловая скорость в положении равновесия согласно уравнению (15): 9

10 ω x = ω, с другой стороны из (17) и (18) получаем следующее выражение для эффективного углового размера спицы Δφ n : ω x = Δφ Δt = ω Δt = Δφ ω Δφ n = Δt () Формула () выполняется точно, если трение отсутствует, при наличии слабого вязкого трения можно сделать поправку первого приближения. В приложении 8. показан вывод с учётом этой поправки выражение (6). 1

11 8 Приложение. 8.1 Учёт вязкого сопротивления при движении маятника. Исходные выражения для связи скорости и угла: φ cos(φ) = β φ dφ + Const cos( ) = Const φ = cos(φ) cos(φ x) β φ dφ φ = 4(sin( ) sin( φ ) ) β φ dφ Обозначим через φ _ скорость маятника без трения (нулевое приближение), и через φ _1 скорость маятника с учётом трения в первом приближении:. φ _ = sin( ) sin( φ ) 1. φ _1 = 4(sin( ) sin( φ ) ) β φ _dφ Далее, для учёта трения, преобразуем интеграл скорости по углу, с учётом нулевого приближения скорости: sin ( ) sin ( φ ) dφ ( ) = sin 1 sin( φ ) dφ (1) sin( ) Максимальный угол отклонения обозначен через. Применим замену переменной: sin( φ ) = sin(θ), dφ = sin( sin( ) 1 sin(θ) dθ ) 1 sin( ) sin(θ) И дальнейшие преобразования интеграла (1) принимают вид: sin( ) sin( π/ ) π/ π/ 1 sin(θ) sin( φ x 1 sin(θ) ) 1 sin( ) sin(θ) dθ = 1 sin(θ) π/ 1 sin( ) sin(θ) dθ = 1 sin( 1 sin( ) sin(θ) dθ+(sin( ) 1) π/ ) sin(θ) + sin( ) 1 1 sin( dθ = ) sin(θ) dθ 1 sin( ) sin(θ) = E( π sin( ) ) + (sin( ) 1) F ( π sin( ) ) Здесь выделены известные специальные математические функции: F ( π sin( ) ) = K(sin( ) ) = π/ dθ — полный эллиптический интеграл первого 1 sin( ) sin(θ) рода 11

12 E( π sin( ) ) = E(sin( ) ) = π/ 1 sin( ) sin(θ) dθ- полный эллиптический интеграл второго рода Существует относительно простой способ численного вычисления этих функций, не требующий явного вычисления интегралов. Используя полученное выражение для интеграла (1) можно записать в первом приближении скорость маятника φ — до и φ 1 — после текущего максимального отклонения: φ = 4 sin( ) + 8β (E(sin( ) ) + (sin( ) 1) K(sin( ) )) φ 1 = 4 sin( ) 8β (E(sin( ) ) + (sin( ) 1) K(sin( ) )) Сумма квадратов скоростей в этом приближении не зависит от параметра затухания: sin( ) = φ + φ 1 () 8 А из разности квадратов скоростей получается оценочное выражение для параметра вязкого трения: φ φ 1 = 16β (E(sin( ) ) + (sin( ) 1) K(sin( ) )) β = φ φ 1 16 (E(sin( ) ) + (sin( ) 1) K(sin( ) )) (3) 8. Полный оборот маятника. Отдельно можно рассмотреть случай, когда максимальная амплитуда отклонения маятника составляет 18 градусов. Тогда интеграл (1) момента трения упрощается и в первом приближении скорость маятника: φ _1 = 4(1 sin( φ ) ) 4β φ π 1 sin( φ ) dφ = 4(1 sin( φ ) ) 8β sin( φ ) 8β = 4(1 sin( φ ))(1+sin(φ )) 8β(1+sin(φ )) = 4(1+sin(φ ))(1 β sin(φ )) А скорость при прохождении положения равновесия: φ 1 = 4 8β cos( φ π )d(φ ) = 4(1 β) (4) — при проходе от амплитуды в 18 градусов до первого пересечения положения равновесия. 1

13 Это выражение можно использовать для определения периода линейной системы: φ = dφ dτ Δφ ω Δt = (1 β) ω = Δφ Δt (1 β) T = 4π Δt (1 β) Δφ (5) А также определения эффективного углового размера спицы: Δφ ω Δt = Δφ n Δt = (1 β) Δφ n = Δφ ω = Δt (1 β) (6) Эффективный угловой размер спицы Δφ n позволяет определять измеряемую угловую скорость сразу в безразмерном виде, а это очень удобно для использования уравнения (11в) для определения амплитуды качаний маятника. 8.3 Численное вычисление полного эллиптического интеграла K(x ) Здесь приводится только алгоритм вычислений, без доказательства. (Обоснование алгоритма: ) На входе в алгоритм подаётся x аргумент функции в диапазоне от до 1 и последовательно выполняются вычисления: 1. a = 1 x. b = 1 a 3. a 1 = a +b 4. b 1 = a b 5. Если a 1 и b 1 различаются меньше чем требуемая точность, то закончить итерации и перейти к п. 7, иначе — перейти к п. 6 6.a = a 1, b = b 1, перейти к п.3 7.K = π 1 x 1 a 1 Вычисленная величина K есть приближённое значение полного эллиптического интеграла K(x ). Итерации этого алгоритма быстро сходятся каждая итерация удваивает количество верных значащих цифр вдвое. При расчёте с помощью калькулятора достаточно провести -3 итерации для получения результата с 4-5 верными знаками после запятой. Ниже показан код С-функции используемой в программе для расчёта значений полных эллиптических интегралов: 13

14 Рис. 5: Код С-функции используемой в программе для расчёта значений полных эллиптических интегралов первого и второго рода 14