Mоделирование процессов и систем в Matlab (966709), страница 20
Текст из файла (страница 20)
п. Создание функций от функций На практике довольно часто возникает необходимость создавать подобные собст- венные процедуры. МАТ?АВ предоставляет такие возможности. Процедура 1еча1 В МАТ1.АВ любая функция (процедура), например с именем Г))й1, может быть выполнена не только с помощью обычного обращения: Ь1.У2...УК1 = ГСНих1.х2,....хп) Ее также можно выполнить посредством специальной процедуры Сееа1: Ь1.у2..ух) = Гека)('Г))И1'.х1.х2....,хп) В этой процедуре имя функции Г))й1 является уже одной вз входных переменных — текстовой (и поэтому помещается между двумя апострофами). Преимуществом второго варианта вызова функции является то, что в этом случае не меняется форма обращения к ней при изменении ее имени, например на ГФ2.
Такой вызов позволяет унифицировать обращение ко всем функциям определенного вида, то есть таким, которые имеют одинаковое число входных и выходных параметров соответствующего типа. Прн этом имя функции (а значит, и сама функция) может быть произвольным и может изменяться при повторных обращениях.
Поскольку при вызове функции с помощью процедуры тека) имя функции рассматривается в качестве одного из входных параметров процедуры, то оно (имя функции) используется как переменная. Поэтому можно оформлять в М-файле обращение к данной функции, не зная еще ее имени. Примеры создания функций от функций Рассмотрим на примерах особенности создания собсгвенных функций от функций. Процедура для метода Рунге-Купа 4-го порядка численного интегрирования ОДУ Допустим, что задана система обыкновенных дифференциальных уравнений (ОДУ) в форме Коши: — = а(у г)) а(у зй где у — вектор переменных состояния системы; г — аргумент (время); 2 — вектор заданных (в общем случае нелинейных) функций, которые, собственно, и определяют конкретную систему ОДУ. Если значение вектора у в момент времени г известно, то общая формула, по которой может быть найден вектор у, значений переменных состояния системы в момент времени г, - г + А (где Ь вЂ” шаг интегрирования)„имеет такой вид: у, -у+АР(уФ).
108 Урок 2 ° Протрамиирование в среде МАТ(АВ Функция Г(у,г), связанная с вектором Е, может принимать разный вид в зависимости от выбранного метода численного интегрирования. Для метода РунгеКутта 4-го порядка выберем следующую ее форму: Р -(й,+зй,+зй,+й,)/В, где )гт = Е(у,с); й, - цу+ Ий,/3, с+ й/3); $(з = Х(у + Нсз Ист/3, г + 2Ь/3); й, -К(у+йй, -йй,-ййпг+й), Создадим М-файл процедуры, осуществляющей эти вычисления; назовем его т((о43: Гипс(топ (свис.уоцг) - гхо43(7ртцп,дд .у) ЙЯКО43 Интегрирование ОДУ нетодон Рунге-Кутта 4-то порядка, 3 правые части которых заданы процедурой 7ртцп $ Входные перененные: $ 7ртып — строка синволов, которая содервит иня процедуры вычисления правых частей ООУ д Обращение: г - тцп(г,у).
где 7р/цп - 'тцп': д à — тенущий нонент вренени; 3 у — вектор текущих значений перененних состояния: д г — вычисленные значения производных г(т) - бу(т)/бг в П вЂ” иат интегрирования; д С вЂ” предыдущий нонент врененн; 3 у — предыдущее значение вектора перененных состояния $ Выходные перененные: д Сонг — новый ионент времени: д уоиг — вычисленное значение вектора у через иат $ Расче~ проиеяуточных значений производных Я1 - теча)(7р(опл,у): 'к2 - /еча)(7р(цп.с + П/З,у + П/3'$т1): хЗ - Гена)(7ртцпД + 2'Д/З.у + пий2 - Ц/3'к1); К4 - Гека)(7р(ипд + п„у + п*(КЗ + Я1 - К2)): $ Расчет новых значений вектора переменных состояния гонг = г + ц; усцг - у + ц*(К1 + Зяй2 + ЗийЗ + й4)/В; з Конец процедуры ЯКО43 Обратите внимание на такие обстоятельства: О обращение к процедуре вычисления правых частей ие конкретизировано; нмя этой процедуры является одной из входных переменных процедуры интегрирования и должно быть конкретизировано лишь при вызове последней; О промежуточные переменные Я являются векторами-строками (как и переменные у и 2, вычисляемые в процедуре правых частей).
Процедура вычисления правых частей Оду маятника Рассмотрим процесс создания процедуры вычисления правых частей ОДУ на примере уравнения движения маятника, точка подвеса которого поступательно перемещается со временем по гармоническому закону: ./ тр+ 47/р+ тттл/(1+ и зщ(тот+ гу ))яптр = — ттфп зш(отг+ в„)сов та Создание функций от функций где: / — момент инерции маятника; >т — коэффициент демпфироваиия; гпф — козффициеит жесткости; и„— амплитуда виброперегрузки точки подвеса маятника в вертикальном направлении; я — амплитуда виброперегрузки в горизонтальном направлении; тр — угол отклонения маятника от вертикали; о> — частота колебаний точки подвеса; е г.
— начальные фазы колебаний точки подвеса в горизоптальиом и вертикальном направлениях. Чтобы составить М-файл процедуры вычисления правых частей задаииой системы ОДУ, прежде всего надо привести исходную систему ОДУ к форме Коши. Для этого введем обозначения: у! >р уз Ф Тогда исходное уравнение движения маятника можно представить в виде совокуп- насти двух дифференциальных уравнений 1-го порядка — = Уз' (у> г>г — = [-тай(яг»» Б>п (о> г + е» )сов (У> ) — гтуг — тпй[[1+ и!»г 5>п (о>г + ед ) зш(у> )) / ). >туг Сравнивая полученную систему с общей формой уравнений Коши, мож>ю сде- лать вывод, что г> =уз; гг = [-тлф п~ з(п(о>1+ е )соз(У> ) — )(Уг — птф[1+ т> г зш(о>г е ед )[ып(У> )) //.
Вычисление именно этих двух функций должно выполняться в процедуре вычислеиия правых частей. Назовем будущую процедуру гИО. Выходной переменной в ией будет вектор г=(г). 22), а входными переменными — момент времени 1 и вектор у-(у1.у2). Определеииую сложность представляет то, что постоянные козффициеиты, находящиеся в правых частях, иельзя передать в процедуру через ее ззголоаок Поэтому обьедииим зти коэффициенты в вектор К=(,),й,од),поу.
пвх,оо.еу,ех) и отнесем этот вектор к категории глобальных: д)оЬа1 К. Тогда М-файл будет иметь следующий вид: гцпсгтоп г гиО(г.у): г Процедура правых частей уравнения физического иаятника г Осунествпяет расчет вектора г" производных г от вектора "у" переиенных состояния по форнулан: г г(1) " у(2): г г(2) - (-ад!*овх'втп(ов 1 + ех)"сов(у(1)) — а"у(2)-..
г ад)*( 1 + пау"зтп(ояяг + еу))"зтп(у(1)))гц г козффициенти передаются в процедуру через гпобапьныд вектор г К - (З.д,зд).оту.пох.ов.еу.ех) О1оЬа! К 7(1) у(2); г(2) " (-К(3)"К(5)*зтп(К(6)ьг + К(З))"соз(у(1)) - К(2)"у(2) †.. К(З)*(1 + К(4)*юп(К(б)яг + К(7)))язтп(у(1)))>К(1); г конец процедуры гиО ыо Урок 2 ° Программирование а среде МАТ(ДВ При использовании этой процедуры следует помнить, что в тексте программы предварительно должен быть объявлен глобальный вектор К с помощью служебного слова д1ОЬа1, а потом определены все восемь его элементов. Данную процедуру можно несколько усложнить, сгруппировав вычисления всех внешних моментов сил, кроме момента сил тяготения, н оформив их как отдельную процедуру.
Для этого сначала преобразуем исходное уравнение, записывая его в таком виде: ~р" + зш~р = Я(т, фч йт'), где штрих обозначает производную по безразмерному времени т = итог, (Я а 5(т,д, ~р') — некоторую заданную функцию безразмерного времени, угла поворота маятника и его безразмерной скорости вар (р = —. гтт В рассматриваемом случае эта функция приобретает такой вид: Я(Кфчяр') = — 2~1р' — (нмя яп(ъ"'с+ ея)сов яр+ Вин яп(от+ еу )51пф), причем безразмерные величины г", и т определяются следующими выражениями: Я го 2 ч(тгтдт' ) ота такая безразмерная форма представления уравнений является предпочтительной, поскольку позволяет сократить количество параметров (в нашем случае вместо трех размерных параметров.)', (т и ац( остался один — т",), а также получить решение уравнения в более общей форме.
Вынесение вычисления моментов внешних действий в отдельную процедуру позволяет также сделать процедуру вычисления правых частей уравнения маятника более общей, если обращение к процедуре вычисления моментов тоже осуществлять через функцию Кеча1. Создадим процедуру Моягв1, которая будет вычислять моменты сил, действующих на маятник: топот)оп и мовтм1(т.у) $ Вичисление клиентов сип. действуощнх на физический иаятник $ Осуществляет расчет нонента "в" сил К по форнуле: $ и - -2"Йяу(2) - (пвхяятп(поят + ех)ясов(у(1)) +... К + пв(тпятп(почт + еу)чзтп(у(1)) К Козффициенти передается в процедуру через глобальний вектор К КМ1 - Иг.пиу.пих.пц.еу.ех) д1оЬа1 КМ1 В - -2ЯКМ1(1)"у(2) — (КМ1(З)чято(КМ1(4)'С + КМН6))ЯСОЗ(у(1)) +...
(падании функций ат функций 111 КМН2)шип(КМ1(4)*С + КМ1(5))вип(у(1))): $ Конец процвдуры МатРМ1 Теперь следует перестроить процедуру правых частей. Назовем этот вариант ГМ1т йасттоп т - ЕМ1(шртмп.т.у) т процедура правих частей уРавнения Физического наятника $ Осуцествляет расчет вектора з" производных д векторов "у" перененных состояния по Ворнулан: т(1) - у(2); т т(2) -втп(у(1)) + 5(т.у) т Входные параметры: т шрттю — иня процедуры 5(т.у); Ф шр(ыт - '5'. Ф С вЂ” текуций ноиент вренени: т у — текуцее значение вентора переиеннмх состояния $ Выхаюае паранетрн: д т — вектор значений производных от перененных состояния т(1) у(2); И2) " -ятп(у(1)) + теча)(шртмл,т.у): т конец процедуры ГМ1 Поскольку вид обращения к процедуре вычисления правых частей изменился (добавлена новая входная переменная — имя процедуры вычисления моментов), необходимо перестроить также процедуру численного метода.
Назовем ее КК043щ: йюсттап (том(.улиц - гко43ш(2ртоп. Мрйю.П.С.у) ВКК043ш Интегрирование ОДу нетодон Рунге-Кутта 4-го порядка. $ правые части котормх заданы проыедурани тртып и Мртып д Входные параметры: й тртмп — строка синволов. которая совершит иня прокедуРы 3 вычисления правых частей Ойу т Обращение: з - йю(мртгю,т,у). д где тртып - 'тта '.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
