Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
h1=text(0.6,0.7,sprintf(' tet0 = %g ',D2(2)),'FontSize',10); h1=text(0.2,0.66,sprintf(' fit0 = %g ',D2(3)),'FontSize',10); h1=text(-0.2,0.62,sprintf('psm = %g ',D2(4)),'FontSize',10); h1=text(0.6,0.62,sprintf(' tem = %g ',D2(5)),'FontSize',10); h1=text(0.2,0.58,sprintf(' fim = %g ',D2(6)),'FontSize',10); h1=text(-0.2,0.54,sprintf('omps = %g ',D2(7)),'FontSize',10); h1=text(0.6,0.54,sprintf(' omte = %g ',D2(8)),'FontSize',10); 2.5. Графическое оформление результатов 117 h1=text(0.2,0.5,sprintf('omfi = %g ',D2(9)),'FontSize',10); h1=text(-0.2,0.46,sprintf('eps = %g ',D2(10)),'FontSize',10); Рис. 2.10 h1=text(0.6,0.46,sprintf(' ete = %g ',D2(11)),'FontSize',10); h1=text(0.2,0.42,sprintf('efi=%g ',D2(12)),'FontSize',10); h1=text(0,0.35,'Интегрирования','FontSize',10,'FontUnderline','on'); h1=text(0,0.3,sprintf('h = %g ',D5(1)),'FontSize',10); h1=text(0,0.25,sprintf('hpr = %g ',D5(2)),'FontSize',10); h1=text(0,0.2,sprintf('t = %g ',D5(3)),'FontSize',10); h1=text(0,0.15,sprintf('tfinal = %g ',D5(4)),'FontSize',10); h1=text(-0.3,0.12,'------------------------------------------------------','FontSize',10); tm=fix(clock); Tv=tm(4:6); h1=text(-0.2,0.08,['Программа ' sprogram],'FontSize',10); h1=text(-0.3,0.04,['Расчеты провел ' sname],'FontSize',10); h1=text(-0.3,0,[sprintf(' %g :',Tv) ' ' date],'FontSize',10); h1=text(-0.3,-0.04,'-----------------------------------------------------','FontSize',10); h1=text(-0.3,-0.08,'Ukraine, KPI, cath. PSON','FontSize',10); Выполнение его приводит к появлению в окне фигуры изображения, пред- ставленного на рис. 2.10. 2.6. Создание функций от функций 118 2.6. Создание функций от функций Некоторые алгоритмы являются общими для функций определенного типа. Поэтому их программная реализация едина для всех функций этого типа, но тре- бует использования алгоритма вычисления конкретной функции. Последний ал- горитм может быть зафиксирован в виде определенной файл-функции. Чтобы первый, более общий алгоритм был приспособлен для любой функции, нужно, чтобы имя этой функции было некоторой переменной, принимающей определен- ное значение (текстового имени файла-функции) только при обращении к основ- ному алгоритму. Такие функции от функций уже рассматривались в разд. 2.1. К ним принад- лежат процедуры: - вычислений интеграла от функции, которые требуют указания имени М- файла, содержащего вычисления значения подынтегральной функции; - численного интегрирования дифференциальных уравнений, использование которых требует указания имени М-файла, в котором вычисляются правые части уравнений в форме Коши; - алгоритмов численного вычисления корней нелинейных алгебраических уравнений (нулей функций), которые нуждаются в указания файла-функции, нуль которого отыскивается; - алгоритмов поиска минимума функции, которую, в свою очередь, надо за- давать соответствующим М-файлом и т.п. На практике довольно часто возникает необходимость создавать собствен- ные процедуры такого типа. MatLAB предоставляет такие возможности. 2.6.1. Процедура feval В MatLAB любая функция (процедура), например, с именем FUN1, может быть выполнена не только с помощью обычного обращения: [y1,y2,... ,yk] = FUN1(x1,x2,... ,xn), а и при помощи специальной процедуры feval: [y1,y2,... ,yk] = feval(‘FUN1',x1,x2,... ,xn), где имя функции FUN1 является уже одной из входных переменных (текстовой - и поэтому помещается в апострофы). Преимуществом вызова функции во второй форме является то, что этот вы- зов не изменяет формы при изменении имени функции, например, на FUN2. Это позволяет унифицировать обращение ко всем функциям определенного вида, т. е. таким, которые имеют одинаковое число входных и выходных параметров опре- деленного типа. При этом имя функции (а значит, и сама функция, которая ис- пользуется) может быть произвольным и изменяться при повторных обращениях. Так как при вызове функции с помощью процедуры feval имя функции рас- сматривается как один из входных параметров процедуры, то его (имя функции) можно использовать как переменную и оформлять в М-файле обращения к ней, не зная еще конкретного имени функции. 2.6. Создание функций от функций 119 2.6.2. Примеры создания процедур от функций Процедура метода Рунге-Кутта 4-го порядка численного интегрирования ОДУ Пусть задана система обыкновенных дифференциальных уравнений (ОДУ) в форме Коши: d dt t y Z y = ( , ) , где y - вектор переменных состояния системы; t - аргумент (время); Z - вектор заданных (в общем случае – нелинейных) функций, которые, собственно, и опре- деляют конкретную систему ОДУ. Если значение вектора y в момент времени t известно, то общая формула, по которой может быть найден вектор yout значений переменных состояния системы в момент времени tout = t + h (где h - шаг интегрирования), имеет вид: yout = y + h*F(y,t). Функция F(y,t) связана с вектором Z и может приобретать разный вид в за- висимости от выбранного метода численного интегрирования. Для метода Рунге- Кутта 4-го порядка выберем такую ее форму: F k k k k 1 2 3 4 = + ⋅ + ⋅ + ( ) 3 3 / 8 / ); 3 , где k Z y 1 = ( , ); t k Z y k 1 2 3 3 = + ⋅ + ( / , h t h k Z y k k 1 3 2 3 2 = + ⋅ − ⋅ + ( / , h h t h / ); k Z y k k k 1 4 3 2 = + ⋅ − ⋅ − ⋅ + ( , h h h t ) h . Создадим М-файл процедуры, которая осуществляет эти вычисления, на- звав его "rko43": function [tout,yout] = rko43(Zpfun,h,t,y) %RKO43 Интегрирование ОДУ методом Рунге-Кутта 4-го порядка, % правые части которых заданы процедурой Zpfun. % Входные переменные: % Zpfun - строка символов, который содержит имя процедуры % вычисления правых частей ОДУ % Обращение: z = fun(t,y), где Zpfun = 'fun' % где t - текущий момент времени % y - вектор текущих значений переменных состояния % z - вычисленные значения производных z(i) = dy(i)/dt. % h - шаг интегрирования % t - предыдущий момент времени % y - предыдущее значение вектора переменных состояния. % Выходные переменные: % tout - новый момент времени % yout - вычисленное значение вектора y через шаг интегрирования % Расчет промежуточных значений производных k1 = feval(Zpfun, t, y); k2 = feval(Zpfun, t+h/3, y+h/3*k1); k3 = feval(Zpfun, t+2*h/3, y+h*k2-h/3*k1); k4 = feval(Zpfun, t+h, y+h*(k3+k1-k2)); 2.6. Создание функций от функций 120 % Расчет новых значений вектора переменных состояния tout = t + h; yout = y + h*(k1 + 3*k2 + 3*k3 + k4)/8; % Конец процедуры RKO43 Обратите внимание на такие обстоятельства: - обращение к процедуре вычисления правых частей не конкретизированно; имя этой процедуры входит в число входных переменных процедуры интегриро- вания и должно быть конкретизировано лишь при вызове последней; - промежуточные переменные k являются векторами-строками (так же, как и переменные 'y' и 'z', вычисляемые в процедуре правых частей). Процедура правых частей ОДУ маятника Рассмотрим процесс создания процедуры вычисления правых частей ОДУ на примере уравнения маятника, точка подвеса которого поступательно переме- щается со временем по гармоничному закону: , 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 ε - начальные фазы колебаний (по уско- рениям) точки подвеса в горизонтальном и вертикальном направлениях. Чтобы составить М-файл процедуры вычисления правых частей заданной системы ОДУ, прежде всего надо привести исходную систему ОДУ к форме Ко- ши. Для этого введем обозначения: y1 = ϕ ; y2 = ϕ & . Тогда исходное уравнение маятника можно представить в виде совокупно- сти двух дифференциальных уравнений 1-го порядка: J y t n mgl y R y t n mgl dt dy y dt dy y my x mx / )} 1 sin( )] sin( 1 [ 2 ) 1 cos( ) sin( { 2 ; 2 1 ⋅ + ⋅ ⋅ + ⋅ − − ⋅ − ⋅ + ⋅ ⋅ ⋅ − = = ε ω ε ω Сравнивая полученную систему с общей формой уравнений Коши, можно сделать вывод, что 2.6. Создание функций от функций 121 J y t n mgl y R y t n mgl z y z y my x mx / )} 1 sin( )] sin( 1 [ 2 ) 1 cos( ) sin( { 2 ; 2 1 ⋅ + ⋅ ⋅ + ⋅ − − ⋅ − ⋅ + ⋅ ⋅ ⋅ − = = ε ω ε ω Именно вычисление этих двух функций и должно происходить в процедуре правых частей. Назовем будущую процедуру fm0. Выходной переменной в ней будет вектор z = [z1 z2], а входными - момент времени t и вектор y = [y1 y2]. Не- которую сложность представляет то, что постоянные коэффициенты в правых частях нельзя передать в процедуру через ее заголовок. Поэтому объединим их в вектор коэффициентов К = [J, R, mgl, nmy, nmx, om, ey, ex] и отнесем этот вектор к категории глобальных global K. Тогда М-файл будет иметь вид: function z = FM0(t,y); % Процедура правых частей уравнения Физического Маятника. % Осуществляет расчет вектора “z” производных % от вектора "y" переменных состояния по формулам: % z(1)=y(2); % z(2)=(-mgl*nmx*sin(om*t+ex)*cos(y(1))-R*y(2)-... % mgl*(1+nmy*sin(om*t+ey))*sin(y(1)))/J, % Коэффициенты передаются в процедуру через глобальный вектор % К=[J,R,mgl,nmy,nmx,om,ey,ex] global K z(1) = y(2); z(2) = (-K(3)*K(5)*sin(K(6)*t+K(8))*cos(y(1)) - K(2)*y(2) -... K(3)*(1+K(4)*sin(K(6)*t+K(7)))*sin(y(1)))/K(1); % Конец процедуры FM0 При использовании этой процедуры следует помнить, что в тексте про- граммы предварительно должен быть объявлен глобальный вектор К с помощью служебного слова global, а потом определены все 8 его элементов. Эту процедуру можно несколько усложнить, группируя вместе вычисление всех внешних моментов сил, кроме момента сил тяготения, и оформляя их как от- дельную процедуру. Для этого сначала превратим начальное уравнение, записы- вая его в виде: ) , , ( sin ϕ ϕ τ ϕ ϕ ′ = + ′′ S , где штрих - обозначение производной по безразмерному времени τ ω = ⋅ 0 t , J mgl = 0 ω , а через ) , , ( ϕ ϕ τ ′ S обозначена некоторая заданная функция безразмерного вре- мени, угла поворота маятника и его безразмерной скорости τ ϕ ϕ d d = ′ В рассматриваемом случае эта функция приобретает такой вид: 2.6. Создание функций от функций 122 ], sin ) sin( cos ) sin( [ 2 ) , , ( ϕ ε τ ν ϕ ε τ ν ϕ ζ ϕ ϕ ⋅ + ⋅ ⋅ + ⋅ + ⋅ ⋅ − − ′ ⋅ ⋅ − = ′ y my x mx n n t S причем безразмерные величины ζ и ν определяются выражениями: J mgl R ⋅ ⋅ = 2 ζ ; 0 ω ω ν = Такая безразмерная форма представления уравнений является предпочти- тельной и удобной, так как позволяет сократить количество параметров (в нашем случае вместо трех размерных параметров J, R и mgl остался один - ζ ), а также представлять решение уравнения в более общей форме. Вынесение вычисления моментов внешних действий в отдельную вычисли- тельную процедуру позволяет также сделать процедуру правых частей уравнения маятника более общей, если обращение к процедуре вычисления моментов осу- ществлять тоже через функцию feval. Создадим процедуру MomFm1, которая будет вычислять моменты сил, ко- торые действуют на маятник: 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 Теперь следует перестроить процедуру правых частей. Назовем этот вари- ант FM1: function z = FM1(mpfun,t,y); % Процедура правых частей уравнения Физического Маятника. % Осуществляет расчет вектора “z” производных % векторов "y" переменных состояния по формулам: % z(1)=y(2); % z(2)= - sin(y(1)) +S(t,y), % входные параметры: % mpfun - имя процедуры S(t,y) % mpfun = 'S'; % t - текущий момент времени; % y - текущее значение вектора переменных состояния; % выходные параметры: % z - вектор значений производных от переменных состояния. z(1) = y(2); z(2) = - sin(y(1)) + feval(mpfun,t,y); % Конец процедуры FM1 2.6. Создание функций от функций 123 Так как вид обращения к процедуре правых частей изменился (добавлена новая входная переменная - имя процедуры вычисления моментов), необходимо также перестроить и процедуру численного метода. Назовем ее RKO43m: function [tout,yout] = rko43m(Zpfun,Mpfun,h,t,y) %RKO43m Интегрирование ОДУ методом Рунге-Кутта 4-го порядка, % правые части которых заданы процедурами Zpfun и Mpfun. % Входные параметры: % Zpfun - строка символов, который содержит имени процедуры % вычисление правых частей ОДУ. % Обращение: z = fun(Mpfun,t,y), % где Zpfun = 'fun', % Mpfun - строка с именем процедуры, к которой % обращается процедура fun; % t - текущий момент времени % y - вектор текущих значений переменных состояния % z - вычисленные значения производных z(i) = dy(i)/dt. % h - шаг интегрирования % t - предшествующий момент времени % y - предшествующее значение вектора переменных состояния. % Выходные параметры: % tout - новый момент времени % yout - новое значение вектора переменных состояния % через шаг интегрирования % Расчет промежуточных значений производных k1 = feval(Zpfun, Mpfun,t, y); k2 = feval(Zpfun, Mpfun, t+h/3, y+h/3*k1); k3 = feval(Zpfun, Mpf un, t+2*h/3, y+h*k2-h/3*k1); k4 = feval(Zpfun, Mpfun, t+h, y+h*(k3+k1-k2)); % Расчет новых значений вектора переменных состояния tout = t + h; yout = y + h*(k1 + 3*k2 + 3*k3 + k4)/8; % Конец процедуры RKO43m Такая форма представления процедуры вычисления правых частей диффе- ренциальных уравнений неудобна. Во-первых, процедуру вида FM1 нельзя ис- пользовать при интегрировании процедурами MatLAB ode23 и ode45 (последние требуют, чтобы в процедуре правых частей было только два входных параметра, а в процедуре FM1 их три). Во-вторых, такая форма вызовет необходимость созда- ния новых М-файлов методов численного интегрирования. Этого можно избегнуть, представив имя дополнительной функции Mpfun на глобальную переменную. Тогда процедура правых частей может быть записана так: |