Картинки из сети, качество желает лучшего, но они достаточно точно отражают суть опыта по визуализации фигур. Зри в корень – основа мудрости поколений.
Немного истории
Ещё в школе на уроках физики я вглядывался в осциллограф, на экране которого, сменяя друг друга, появлялись разные фигуры: сначала простые – линия, парабола, круг, эллипс, потом фигуры становились всё более насыщенные непрерывными волнообразными линиями, напоминающие мне кружева. Автором этого кружевного дива был Жюль Антуан Лиссажу французский физик, член — корреспондент Парижской АН (1879) [1]. Сами фигуры — это замкнутые траектории, прочерчиваемые точкой, совершающей одновременно два гармонических колебания в двух взаимно перпендикулярных направлениях [2]. Думаю, что в те далёкие от современности годы основной заслугой Жюля, кроме конечно накопленных опытом знаний математики и физики, была простая механическая визуализация этих фигур подручными средствами. Захотелось конструировать подобно Жулю максимально просто и наглядно, реализовать его идеи применительно к современной задаче линейных измерений. Но сделать это путём математического моделирования с графической визуализацией его результатов на Python. Но сначала рассмотрим классический вариант [3] построения фигур.
Какими должны быть фигуры Лиссажу
Для этого воспользуемся системой уравнений, описывающих фигуры:
x(t), y(t) в общем случае зависящие от времени гармонические колебания вдоль взаимно перпендикулярных плоскостей, частоты b, a и начальная фаза d. Для анализа фигур в вычислениях принимают постоянным модуль разности частот |b — a| = 1. Будем рассматривать отношение круговых частот b / a и начальную фазу d. Имеем для линии A = B d = 0, окружности 
, и параболы 
. Основные отношения частот, удовлетворяющие условию, занесём во вложенный список m=[[0],[2,2],[2,1],[1,2],[3,2],[3,4],[5,4],[5,6],[9,8]].
#!/usr/bin/env python #coding=utf8 import numpy as np from numpy import sin,pi import matplotlib.pyplot as plt m=[[0],[2,2],[2,1],[1,2],[3,2],[3,4],[5,4],[5,6],[9,8]]# отношение круговых частот for i in m: if i[0]==0: a=1 x=[sin(a*t) for t in np.arange(0.,2*pi,0.01)] y=[sin(a*t) for t in np.arange(0.,2*pi,0.01)] plt.plot(x, y, 'r')# график для линии plt.grid(True) plt.show() else: a=i[0] b=i[1] d=0.5*pi x=[sin(a*t+d) for t in np.arange(0.,2*pi,0.01)] y=[sin(b*t) for t in np.arange(0.,2*pi,0.01)] plt.plot(x, y, 'r') # график для различных отношений a/b #круговых частот plt.grid (True) plt.show() Результат не привожу, отдельные фигуры не впечатляют. Хочу коллаж из «кружев».
#!/usr/bin/env python #coding=utf8 import numpy as np from numpy import sin,pi import matplotlib.pyplot as plt m=[[3,4],[5,4],[5,6],[9,8]] # отношение круговых частот plt.figure(1) for i in m: a=i[0] b=i[1] d=0.5*pi x=[sin(a*t+d) for t in np.arange(0.,2*pi,0.01)] y=[sin(b*t) for t in np.arange(0.,2*pi,0.01)] if m.index(i)==0: plt.subplot(221) plt.plot(x, y, 'k') # график для различных отношений a/b круговых частот plt.grid(True) elif m.index(i)==1: plt.subplot(222) plt.plot(x, y, 'g') plt.grid(True) elif m.index(i)==2: plt.subplot(223) plt.plot(x, y, 'b') plt.grid(True) else: plt.subplot(224) plt.plot(x, y, 'r') plt.grid(True ) plt.show() И вот они «кружева».
Что нельзя отнести к фигурам Лиссажу по определению о их замкнутости
Зачем нам |b — a| = 1, “за флажки!” попробуем например так m=[[1,3],[1,5],[1,7],[1,9]]
На втором графике при m=0,2 получена незамкнутая траектория, которая по определению не является фигурой Лbссажу.
В поисках механических аналогов
Поищем аналогии фигур в измерительной технике и вот вибрационный уровнемер с резонатором в виде эллиптической трубки [4].
Упруго закреплённая трубка эллиптического сечения с помощью систем возбуждения 5,6,7 совершает автоколебания в одной плоскости, а с помощью систем 8, 9, 10 в другой плоскости перпендикулярной первой. Трубка колеблется в двух взаимно перпендикулярных плоскостях с разными частотами близкими к собственным. Масса трубки зависит от уровня заполняющей её жидкости. С изменением массы меняются и частоты колебаний трубки, которые и являются выходными сигналами уровнемера. Частоты несут дополнительную информацию о мультипликативных и аддитивных дополнительных погрешностях, компенсируемых при обработке частот микропроцессором 11.
Условия адекватного моделирования
Для более-менее корректной привязки фигур Лиссажу к работе упомянутого уровнемера, следует учесть следующие обстоятельства. Во-первых, закреплённая одним концом трубка эллиптического сечения — это колебательная система с распределёнными параметрами, что сильно усложняет анализ её колебаний. Во-вторых, отношение частот колебаний трубки не может изменяться произвольно, оно зависит от эллипсности сечения и допустимых зазоров в системе возбуждения колебаний. Для отношения частот можно получить простое соотношение.
К чему принадлежат переменные, a, b, a0, b0 ясно из рисунка и кроме того формула для циклической частоты осциллятора известна из школьного курса физики. Для «реализации на Python в последнее отношение введём толщину стенки и показатель эллипсности внутреннего сечения трубки, тогда вместо четырёх переменных получим три.
#!/usr/bin/env python #coding=utf8 import numpy as np from numpy import sqrt import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.family'] = 'fantasy' mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial' d=0.5 a=9 x=[w for w in np.linspace(0.8,0.95,15)] y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in np.linspace(0.8,0.95,15)] plt.plot(x, y, 'r', label='Толщина стенки трубки в мм. -- %s' %str(d)) d=0.7 y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in np.linspace(0.8,0.95,15)] plt.plot(x, y, 'b',label='Толщина стенки трубки в мм.-- %s' %str(d)) d=1.0 y=[sqrt(x**-2*((1+d/a)**3*(1+d/(x*a))-1)/((1+d/a)*(1+d/(x*a))**3-1)) for x in np.linspace(0.8,0.95,15)] plt.plot(x, y, 'g', label='Толщина стенки трубки в мм.-- %s' %str(d)) plt.ylabel('Отношение частот колебаний эллиптической трубки') plt.xlabel('Отношение длин малой и большой полуосей') plt.title('Определение допустимого диапазона для отношения частот') plt.legend(loc='best') plt.grid(True) plt.show() В результате работы программы получим график.
График построен для малой внутренней полуоси в 9 мм. Для конструктивно допустимого отношения малой к большой полуоси сечения в диапазоне от 0.8 до 0.95. Это основной фактор влияния на отношение частот, которое изменяется от 1.18 до 1.04. Толщина стенки влияет незначительно. Теперь у нас есть диапазон отношений и ним можно воспользоваться для дальнейшего моделирования.
Формы колебаний вертикальной оси трубки
Что касается распределённых механических параметров консольной трубки, то они при помощи равенства собственных частот и импеданса могут быть приведены к сосредоточенной массе жёсткости и демпфированию. Кроме того, для определения форм изгибных колебаний консольной трубки можно получить выражение для распределённых параметров. Уравнение для форм – балочные функции имеет вид:
где 
— корни уравнения:
Следует отметить что, не смотря на большое количество публикаций о формах и частотах колебаний консольного стержня, балки или трубки уравнения (4) нигде не приводяться, только рисунки без координат. Поэтому уравнение (4), я вывел через условия на концах и балочные функции, проверил по корням (5) и расположению узлов. Однако это тривиальное уравнение, о котором просто забыли.
#!/usr/bin/env python #coding=utf8 from scipy.optimize import * import numpy as np from numpy import pi,cos,cosh,sin,sinh import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.family'] = 'fantasy' mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial' d=[] for i in range(0,4): x=brentq(lambda x:cosh(x)*cos(x)+1,0+pi*i,pi+pi*i) p=round(x,3) if p not in d: d.append(p) x=[w for w in np.linspace(0,1,100)] k=d[0] z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x) )for x in np.linspace(0,1,100)] plt.plot(z, x, 'g', label='Первая форма для корня - %s' %str(k)) k=d[1] z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x)) for x in np.linspace(0,1,100)] plt.plot(z, x, 'b', label='Вторая форма для корня - %s' %str(k)) k=d[2] z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x)) for x in np.linspace(0,1,100)] plt.plot(z, x, 'r', label='Третья форма для корня - %s' %str(k)) plt.title('Первые три формы изгибных колебаний осевой линии трубки') plt.xlabel(' Координата вдоль оси OX ') plt.ylabel(' Координата положения осевой линии трубки вдоль оси OZ ') plt.legend(loc='best') plt.grid(True) plt.show() В результате работы программы получим график построенный с учётом вертикального положения трубки.
На графике координата осевой линии приведена к длине трубки, а амплитуда нормирована. Положение узлов колебаний трубки относительно места её крепления в точности соответствует теории колебаний.
По каким траекториям движется конец трубки
Последнее препятствие — сложность получения осмысленного численного решения дифференциальных уравнений колебаний, при условии варьирования несколькими параметрами одновременно. Тут на помощь пришли две мои статьи о колебательном звене на Python [5,6], в которых приведена методика получения точных символьных решений дифференциальных уравнений.
Запишем два условно независимых уравнения для колебаний трубки в плоскости OX и OY с разными частотами a и b отношение между которыми выбрано из ранее установленного диапазона. Остальные параметры выбраны во правильной взаимосвязи, но произвольно для лучшей демонстрации результата.
Здесь введены следующие обозначения (для упрощения без индексов).
─ приведенная амплитуда силы, 
─ коэффициент затухания, 
─ собственная частота колебаний системы, m ─ сосредоточенная масса одинаковая для обоих уравнений, 
─ сосредоточенные коэффициенты демпфирования, разные из-за разных амплитуд, а следовательно разных зазорах в системах возбуждения колебаний, 
─ разные жёсткости из-за эллиптичности сечения трубки.
import numpy as np from sympy import * from IPython.display import * import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.family'] = 'fantasy' mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Arial' def solution(w,v,i,n1,n2,B,f,N): t=Symbol('t') var('t C1 C2') u = Function("u")(t) de = Eq(u.diff(t, t) +2*B*u.diff(t) +w**2* u, f*sin(w*t+v)) des = dsolve(de,u) eq1=des.rhs.subs(t,0) eq2=des.rhs.diff(t).subs(t,0) seq=solve([eq1,eq2],C1,C2) rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2])]) g= lambdify(t, rez, "numpy") t= np.linspace(n1,n2,N) plt.figure(1) if i==1: plt.subplot(221) plt.plot(t,g(t),color='b', linewidth=3,label='x=%s*sin(%s*t+%s)' %(str(f),str(w),str(v))) plt.legend(loc='best') plt.grid(True) else: plt.subplot(222) plt.plot(t,g(t),color='g', linewidth=3,label='y=%s*sin(%s*t+%s)' %(str(f),str(w),str(v))) plt.legend(loc='best') plt.plot(t,g(t),color='r', linewidth=3) plt.grid(True) return g(t) N=1000#Число точек оцифровки временного интервала B=0.2#Установка демпфирования f=1#Установка амплитуды n1=0#Нижняя граница временной развертки n2=20#Верхняя граница временной развёртки w1=5.0#Частота колебаний трубки вдоль оси ОХ w2=10.0#Частота колебаний трубки вдоль оси ОУ v1=0#Начальная фаза при колебании вдоль оси ОХ v2=0#Начальная фаза при колебании вдоль оси ОУ g1=solution(w1,v1,1,n1,n2,B,f,N) g2=solution(w2,v2,2,n1,n2,B,f,N) plt.subplot(223) plt.plot(g1,g2,color='b', linewidth=3,label='w1/w2=%s'%str(w1/w2)) plt.legend(loc='best') plt.grid(True) plt.subplot(224) x=[w for w in np.linspace(0,1,100)] k=1.875 z=[sin(k*x)-sinh(k*x)+((cosh(k) -cos(k))/(sin(k)-sinh(k)))*(cos(k*x)-cosh(k *x) )for x in np.linspace(0,1,100)] plt.plot(z, x, 'g',label='Форма -%s'%str(k)) plt.legend(loc='best') plt.grid(True) plt.show() Программа позволяет менять все параметры модели, например, для:
N=1000, B=0.2, f=1, n1=0, n2=20, w1=5.0, w2=10.0, v1=0, v2=0
Для отношения частот 0.5 переходной процесс множит фигуры. Поставим “ворота” времени n15=0, n2=20, получим.
Снимем” ворота” и введём начальную фазу v2=-pi/2, получим:
С учётом изложенного выше, графики комментарий не требую.
Для интриги
Если эта статья найдёт своих читателей или читатели её найдут, не устрашившись теней прошлого, то я опубликую трёхмерные анимационные графики сложных пространственных колебаний трубки при изменении в ней уровня заполняющей жидкости.
Вместо выводов
Изобретение Жюля Антуана Лиссажу продолжает свой путь во времени, но уже и на Python. Надеюсь, что представленная интерпретация, конечно далёкая от совершенства, позволит продолжить знакомство с работами гениального математика Лиссажу.
Ссылки
- Биографии учёных физиков.
- Что такое фигуры Лиссажу?
- Фигуры Лиссажу.
- Вибрационный уровнемер.А.С.№777455
- Модель колебательного звена с применением символьного и численного решений дифференциального уравнения на SymPy и NumPy.
- Модель колебательного звена в режиме резонансных колебаний на Python.
ссылка на оригинал статьи https://habrahabr.ru/post/327768/
Добавить комментарий