Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
2.6. Создание функций от функций 130 где - массы гравитирующих материальных тел; - радиу- сы-векторы материальных точек относительно их общего центра масс; m m m 1 2 , , 3 3 R R R 1 2 , , 3 1 13 2 3 32 1 2 21 R R R R R R R R R − = − = − = ; ; - радиусы-векторы точек друг относи- тельно друга; - текущие расстояния между точками; R R R 21 32 13 , , γ - гравитаци- онная постоянная. Задание 2.24. Составьте процедуру правых частей дифференциального уравнения лампового генератора: q q q q & && ) ( 2 2 λ χ ω − = + Это так называемое уравнение Ван-дер-Поля. Задание 2.25. Составьте процедуру правых частей дифференциального уравнения движения маятника на подвижном основании с учетом момента сил су- хого трения вдоль оси маятника: ), sin( sin 0 t M M M mgl R J m тр ω ϕ ψ ϕ + + = + + & && где ) (t ϑ ϕ ψ − = , ) (t ϑ - текущий угол поворота основания вокруг оси маятника; - постоянная составляющая момента внешних сил, которые действуют на ма- ятник; - амплитуда гармонической составляющей момента внешних сил; 0 M m M ω - частота изменения момента сил; - момент сил сухого трения по оси маятника. тр M Считать ), ( ψ & sign M M тр − = где ⎪ ⎩ ⎪ ⎨ ⎧ < − = > + = , 0 п , 1 , 0 п , 0 , 0 п , 1 ) ( x ри x ри x ри x sign а также ) sin( ) ( 0 0 ε ω ϑ ϑ ϑ ϑ + + + = t t t m & Задание 2.26. Составьте процедуру отыскания точных частных решений системы линейных неоднородных дифференциальных уравнений по заданной матрице системы ДУ в форме Коши и матрице коэффициентов входных воздействий: A B d dt y A y B x = + , где вектор входных воздействий принять в виде x x B B B o S C = + + sin( ) cos( ). ω ω t t 2.7. Пример создания сложной программы Ранее были рассмотрены основные препятствия, стоящие на пути создания сложных программ, и средства их преодоления. Теперь, учитывая это, попробуем составить и испытать в работе одну из довольно сложных комплексных программ. 2.7. Приклад складної програми 131 2.7.1. Программа моделирования движения маятника Постановка задачи. Пусть требуется создать программу, которая позволи- ла бы моделировать движение физического маятника с вибрирующей точкой под- веса путем численного интегрирования дифференциального уравнения этого дви- жения. Дифференциальное уравнение движения маятника для этой задачи можно принять таким: , cos ) sin( sin )) sin( 1 ( ϕ ε ω ϕ ε ω ϕ ϕ ⋅ + ⋅ ⋅ ⋅ − = = ⋅ + ⋅ ⋅ + ⋅ + ⋅ + ⋅ x mx y my t n mgl t n mgl R J & && где: J - момент инерции маятника; R - коэффициент демпфирования; mgl - опор- ный маятниковый момент маятника; - амплитуда виброперегрузки точки под- веса маятника в вертикальном направлении; my n mx n - амплитуда виброперегрузки в горизонтальном направлении; ϕ - угол отклонения маятника от вертикали; ω - частота колебаний точки подвеса; ε x , ε y - начальные фазы колебаний (по ускоре- ниям) точки подвеса в горизонтальном и вертикальном направлениях. Нужно создать такую программу, которая позволяла бы вычислять закон изменения угла отклонения маятника от вертикали во времени при произвольных, устанавливаемых пользователем значениях всех вышеуказанных параметров ма- ятника и поступательного движения основания, а также при произвольных на- чальных условиях. Вычисления будем осуществлять путем численного интегри- рования с помощью стандартной процедуры ode45 Преобразование уравнения. Для подготовки дифференциальных уравне- ний к численному интегрированию прежде всего необходимо привести эти урав- нения к нормальной форме Коши . Желательно также представить их в безразмер- ной форме. Для представленного уравнения это было сделано в предідущем раз- деле. Запись М-файла процедуры правых частей . Следующим шагом подго- товки программы является написание и запись на диск текста процедуры вычис- ления правых частей полученной системы ДР в форме Коши. Эта процедура была создана в предшествующем разделе и в окончательном варианте она имеет вид: Файл FM2 . m function z = FM2(t,y); % Процедура правых частей уравнения Физического Маятника. % Осуществляет расчет вектора "z" производных вектора % "y" переменных состояния по формулам: % z(1)=y(2); % z(2)=-sin(y(1)) +S(t,y), % Входные параметры: % t - текущий момент времени; 2.7. Приклад складної програми 132 % y - текущее значение вектора переменных состояния; % MPFUN - имя процедуры S(t,y) - глобальная переменная % MPFUN = 'S'; % Выходные параметры: % z - вектор значений производных от переменных состояния global MPFUN z(1) = y(2); z(2) = - sin(y(1))+ feval(MPFUN,t,y); z =z'; % Конец процедуры FM2 Примечание. Обратите внимание на незначительное, но существенное от- личие приведенной процедуры от аналогичной процедуры предыдущего раздела - наличие в конце процедуры операции транспонирования вектора производных. Это обусловлено тем, что процедура ode45 требует, чтобы вектор производных был обязательно столбцом. В качестве дополнительной процедуру, используемой в процедуре FM2, выберем ранее созданную процедуру MomFM1, которую запишем в файл MomFM1. Файл MomFM1 . m function m = MomFM1(t,y); % Вычисление Моментов Сил, которые действуют на Физический Маятник. % Осуществляет расчет момента "m" сил % по формуле: % m =-2*dz*y(2) - (nmx*sin(nu*t+ex)*cos(y(1)) +... % + nmy*sin(nu*t+ey)*sin(y(1)), % Коэффициенты передаются в процедуру через % глобальный вектор КM1=[dz,nmy,nmx,nu,ey,ex] global KM1 m = -2*KM1(1)*y(2)- KM1(3)*sin(KM1(4)*t+KM1(6))*cos(y(1)) -... KM1(2)*sin(KM1(4)*t+KM1(5))*sin(y(1)); % Конец процедуры MomFM1 Очевидно, что в вызывающем Script-файле надо предусмотреть объявление имени дополнительного файла MomFM1 как глобальной переменной MPFUN, а также обеспечить объявление глобальной переменной по имени KM1 и задание значений этого числового массива из пяти элементов. Создание управляющего (главного) Script-файла. Главный файл созда- дим соответственно рекомендациям деления. 2.4.5 : Файл FizMayatn2 . m % FizMayatn2 % Управляющая программа исследования движения физического маятника, % установленного на поступательно вибрирующем основании FizMayatn2_Zastavka k = menu(' Что делать ? ',' Продолжить работу ', ' Закончить работу '); if k==1, while k==1 FizMayatn2_Menu FizMayatn2_Yadro k =menu(' Что делать ? ',' Продолжить работу ', ' Закончить работу '); end 2.7. Приклад складної програми 133 end clear global clear % Конец FizMayatn2 Как видим, программа вызовет три дополнительных Script-файла - FizMa- yatn2_Zastavka, FizMayatn2_Menu и FizMayatn2_Yadro. Поэтому нужно создать еще эти три М-файла. Создание Script-файла заставки. Как отмечалось, этот файл должен со- держать операторы вывода на экран информации об основных особенностях ма- тематической модели, реализованной в программе, и ввода исходных значений параметров этой модели. Ниже приведен текст М-файла FizMayatn2_Zastavka. Файл FizMayatn2_Zastavka . m % FizMayatn2_Zastavka % Часть (вывод заставки на экран) программы FizMayatn2 % Ввод "вшитых" значений sprogram = 'FizMaytn2.м'; sname ='Лазарев Ю.Ф.'; KM1 = [0 0 0 0 0 0]; MPFUN = 'MomFm1'; global KM1 MPFUN tfinal =2*pi*5; fi0 =pi/180; fit0 = 0; clc disp( [' Это программа, осуществляющая интегрирование уравнения ';... ' Физического Маятника при поступательной вибрации точки подвеса ';... ' в форме ';... ' fi" + sin(fi) = - 2*dz*fi'' - ';... ' - nmy*sin(nu*t+ey)*sin(fi) -nmx*sin(nu*t+ex)*cos(fi)';... 'где fi - угол отклонения маятника от вертикали, ';... ' dz - относительный коэффициент затухания, ';... ' nu - относительная частота вибрации точки подвеса, ';... ' nmy, nmx - амплитуды виброперегрузки в вертикальном ';... ' и горизонтальном направлениях соответственно, ';... ' еy,еx - начальные фазы колебаний в вертикальном ';... ' и горизонтальном направлениях соответственно, ';... ' KM1 = [dz,nmy,nmx,nu,ey,ex] - матрица коэффициентов ']) % Конец FizMayatn2_Zastavka В нем осуществляется присваивание исходных ("вшитых") значений всем параметрам заданного дифференциального уравнения, а также параметрам чис- ленного интегрирования - начальным условиям движения маятника и длительно- сти процесса интегрирования. Часть этих параметров объединяется в единый гло- бальный вектор КМ1. Одновременно переменной MPFUN, которая будет использоваться при интегрировании, присваивается значение 'MomFm1'. Создание файла меню. Содержимое файла меню FizMayatn2_Menu приве- дено ниже. Файл FizMayatn2_Menu . m % FizMayatn2_Menu % Часть (осуществляющая диалоговое изменение данных) 2.7. Приклад складної програми 134 % программы FizMayatn2 k=1; while k<10 disp(' ') disp(' Сейчас установлено ') disp([sprintf(' Начальный угол (градусы) = %g', fi0*180/pi),... sprintf(' Начальная скорость = %g', fit0)]) disp(sprintf(' Число периодов = %g', tfinal/2/pi)) KM1 % КМ1=[dz,nmy,nmx,nu,ey,ex] k = menu( ' Что изменять ? ', ... sprintf(' Относительный к-нт затухания = %g', KM1(1)),... sprintf(' Перегрузка (вертикаль) = %g', KM1(2)),... sprintf(' Перегрузка (горизонталь) = %g', KM1(3)),... sprintf(' Относительная частота = %g', KM1(4)),... sprintf(' Фаза (вертикаль) = %g', KM1(5)),... sprintf(' Фаза (горизонталь) = %g', KM1(6)),... sprintf(' Начальный угол (градусы) = %g', fi0*180/pi),... sprintf(' Начальная скорость = %g', fit0),... sprintf(' Количество периодов = %g', tfinal/2/pi),... ' Ничего не изменять '); disp(' ') if k<7, KM1(k) = input( ['Сейчас KM1(',num2str(k),sprintf(') = %g', KM1(k)),... ' Введите новое значение = ']); elseif k==7, fi0 = input([sprintf('Сейчас fi0 = %g градусов', fi0*180/pi),... ' Введите новое значение = ']); fi0 = fi0*pi/180; elseif k==8, fit0 = input([sprintf('Сейчас fit0 = %g', fit0),... ' Введите новое значение = ']); elseif k==9,tfinal=input([sprintf('Сейчас количество периодов = %g', tfinal/2/pi),... ' Введите новое значение = ']); tfinal = tfinal*2*pi; end end % FizMayatn2_Menu Файл осуществляет организацию диалоговой ввода-изменения значений па- раметров физического маятника, движения основания и параметров численного интегрирования в соответствия со схемой, описанной в разд. 2.4.4. Создание файла ядра программы. Основные действия по организации процесса численного интегрирования и выведению графиков сосредоточены в файле FizMayatn2_Yadro: Файл FizMayatn2_Yadro . m % FizMayatn2_Yadro % Часть (осуществляющая основные вычисления) % программы FizMayatn2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1. Подготовка начальных условий %---------------------------------- t = 0; tf = tfinal; y0 =[fi0 fit0]; options = odeset('RelTol',1e-8,'AbsTol',[1e-10 1e-10]); %---------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 2. Организация цикла интегрирование 2.7. Приклад складної програми 135 %---------------------------------- [t,y] = ode45('FM2',[0 tf],y0,options); %---------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 3. Вывод графиков subplot(2,1,2); plot(t/2/pi,y(:,1)*180/pi);grid; title('Отклонение от вертикали','FontSize',14); xlabel('Время (в периодах маленьких собственных колебаний)','FontSize',12); ylabel('Угол в градусах','FontSize',12); subplot(2,4,1:2); plot(y(:,1)*180/pi,y(:,2));grid; title('Фазовая траектория','FontSize',14); xlabel('Угол в градусах','FontSize',12); ylabel('Скорость','FontSize',12); %______________________ % Вывод текстовой информации в графическое окно subplot(2,4,3:4); axis('off'); h1=text(0,1.1,'Моделирование движения физического маятника', 'FontSize', 14, 'FontWeight', 'Bold'); h1=text(0.4, 1,'за уравнением','FontSize',12); h1=text(0,0.9,'fi" + 2*dz*fi'' + [1+nmy*sin(nu*t+ey)]*sin(fi) =','FontSize',14); h1=text(0.55,0.8,' = - nmx*sin(nu*t+ex)*cos(fi)','FontSize',14); h1=text(0,0.7,'за таких значений параметров:','FontSize',12); h1=text(0.45,0.6,sprintf('dz = %g',KM1(1)),'FontSize',12); h1=text(0,0.5,sprintf('nmy = %g',KM1(2)),'FontSize',12); h1=text(0.7,0.5,sprintf('nmx = %g',KM1(3)),'FontSize',12); h1=text(0,0.4,sprintf('ey = %g',KM1(5)),'FontSize',12); h1=text(0.7,0.4,sprintf('ex = %g',KM1(6)),'FontSize',12); h1=text(0.45,0.3,sprintf('nu = %g',KM1(4)),'FontSize',12); h1=text(0,0.2,'и начальных понимал:','FontSize',12); h1=text(0,0.1,[sprintf('fi(0) = %g',fi0*180/pi),' градусов'],'FontSize',12); h1=text(0.7,0.1,sprintf('fi''(0) = %g',fit0),'FontSize',12); h1=text(0,0.05,);------------------------------------------------------------------------------'); h1=text(0,-0.2,);------------------------------------------------------------------------------'); h1=text(-0.05,-0.05,['Программа ',sprogram]); h1=text(0.55,-0.05,'Автор - Лазарев Ю.Ф., каф. ПСОН'); h1=text(0,-0.15,['Выполнил ',sname]); tm=fix(clock); Tv=tm(4:5); h1=text(0.65,-0. 15,[sprintf(' %g:',Tv),' ',date]); % Конец файла FizMayatn2_Yadro Как видим, основные операции включают три главные группы - ввод на- чальных условий, организацию цикла интегрирования и организацию оформления графического окна вывода. Отладка программы. Отладка программы состоит из запуска главного М- файла FizMayatn2, проверки правильности функционирования всех частей про- граммы, внесение корректив в тексты используемых М-файлов до тех пор, пока все запрограммированные действия не будут удовлетворять заданным требовани- ям. Сюда же входят и действия по проверке "адекватности", т. е. соответствия по- лучаемых программой результатов отдельным априорно известным случаям пове- дения исследуемой системы. Очевидно, для такой проверки нужно подобрать не- 2.7. Приклад складної програми 136 сколько совокупностей значений параметров системы, при которых поведение системы является известным из предыдущих теоретических или эксперименталь- ных исследований. Если полученные программой результаты полностью согласу- ются с известными, программа считается адекватной принятой математической модели. В приведенном тексте программы "вшитые" начальные значения парамет- ров отвечают свободному движению маятника без влияния трения. При таких ус- ловиях движение маятника представляет собой незатухающие колебания относи- тельно положения вертикали. Поэтому, если программа работает верно, на графи- ках должны наблюдаться именно такие колебания маятника. Результат работы созданной программы при этих условиях представлен на рис. 2.11. Как видно, в этом отношении программа является адекватной принятой математической моде- ли. Проведение исследований. Созданная программа теперь может быть ис- пользована для моделирования и исследования разнообразных нелинейных эф- фектов, которые наблюдаются у физического маятника при поступательной виб- рации точки его подвеса. На рис. 2.12 - 2.17 продемонстрированные некоторые возможности созданной программы. На рис. 2.12 приведены параметрические колебания маятника, которые мо- гут возникать при вибрации точки подвеса в вертикальном направлении. Из ри- сунка видно, что в этом случае колебания маятника относительно вертикали сна- чала увеличиваются по амплитуде, а потом устанавливаются, причем частота по- стоянных колебаний вдвое меньше частоты вибрации основания и составляет примерно 1,15. Выпрямительный эффект маятника проиллюстрирован на рис. 2.13. В этом случае одновременная вибрация основания в вертикальном и горизонтальном на- правлениях приводит к отклонению среднего положения маятника от вертикали на угол около -5 градусов. Рис. 2.14 иллюстрирует стационарные колебания маятника относительно верхнего положения равновесия, которые могут наблюдаться при интенсивной вертикальной вибрации. Эти колебания при наличии трения затухают, как показан на рис. 2.15, и маятник "застывает" в верхнем положении. |