Решение краевых задач. Лабораторная работа №2б (1253748), страница 2
Текст из файла (страница 2)
значению вектора t.Полная форма обращения к функции ode имеет следующий вид:[tout, xout, varargout] = ode***(fun, tspan, x0, options, varargin).Здесь: fun - имя m-файла, являющегося функцией MATLAB от t и x, в котором вычисляется_____________________________________________________________________________ 4Деменков Н.П. Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________вектор функция f(x,t), т.е. правые части системы дифференциальных уравнений; tspan –вектор, задающий интервал интегрирования [t0 tfinal], t0 – начальное значениеинтервала, tfinal – конечное; могут быть заданы и промежуточные контрольные значениянезависимой переменной, тогда tspan=[t0 t1…tfinal]; x0 – вектор начальных условий (скалярили вектор столбец); tout – вектор-столбец контрольных значений независимой переменной;если используется минимальный вариант tspan, выдаются все значения, которые получались впроцессе численного интегрирования; если tspan содержит и другие значения, кроме t0 и tfinal,то tout=tspan; xout – решения, представленные массивом, в котором каждая строкасоответствует одному элементу в столбце tout; options – аргумент, позволяющийзадавать управляющие параметры; varargin – дополнительные аргументы для вычисленияfun; varargout – дополнительные результаты, возникающие при некоторых вариантахзадания options.Аргумент options позволяет отразить очень много управляющих параметров,которые задаются путем обращения к функции odeset, аргументом которой являетсяпоследовательность пар вида <’параметр’, значение>, в которой название любого изпараметров должно быть заключено в апострофы.Это может быть строка параметров, определяющих значения допустимойотносительной (параметр RelTol) и абсолютной (параметр AbsTol) погрешностейинтегрирования.
Эти параметры можно не указывать, если пользователя устраиваютзначения погрешностей, заданных по умолчанию (1e-3) и (1e-6) соответственно. В противномслучае, перед обращением к процедуре ode следует указать значения погрешностей припомощи процедуры odeset:options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]).Параметр InitialStep аргумента options задает начальный шаг, который поумолчанию выбирается автоматически.Параметр Mass определяет квадратную матрицу масс M, которая являетсямножителем при производной в левой части дифференциального уравненияM(t,x) x =f(t,x).Значением этого параметра является сама матрица, если она заполненаконстантами, или указание на m-файл, содержащий функцию, которая вычисляет матрицумасс, как функцию переменных t и x, например,options=odeset('Mass',@mass).Параметр MassSingular позволяет сообщить, является ли (может ли быть)матрица Mass вырожденной.Параметр Events задает события, наступление которых влияет на ход вычислений.Значением этого параметра является указание на m-файл, содержащий функцию, котораяотслеживает эти события, например,options=odeset('Events',@events).Пример условия окончания процесса интегрированияfunction(value, isterminal, direction)=events (t,x)фиксация времени, когда x(1) проходит через ноль в направлении убывания ипрекращения интегрированияvalue= x(1); - для слежения за обращением в ноль величины x(1); isterminal=1; прекратить интегрирование при value= 0; direction=-1; - при условии убывания value.В этой функции отслеживается обращение в ноль переменной value (необходимымусловием останова является обращение в ноль некоторого выражения).Переменная isterminal задается равной 1, если при достижении value=0интегрирование должно прекращаться, и равной 0, если интегрирование должнопродолжаться.Переменная direction должна принимать одно из трех значений: 1, если valueвозрастает при переходе через ноль; -1, если value убывает при переходе через ноль; 0 прилюбом обращении value в ноль.Все три выходные переменные могут быть векторами одинаковой длины.
В этом_____________________________________________________________________________ 5Деменков Н.П. Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________случае их компоненты соответствуют друг другу: при value(i)=0 интегрированиепрекращается, если isterminal(i)=1 и продолжается, если isterminal(i)=0.
Аналогичноиспользуются и компоненты direction(i).Параметр Jacobian задает в аналитическом виде якобиан f , используемый вxнеявных методах. Значением этого параметра является сама матрица Якоби, если она всязаполнена константами, или указание на m-файл, содержащий функцию, которая вычисляетэту матрицу, как функцию переменных t и x, например,options=odeset('Jacobian',@FJac).Параметр Jacobian имеет смысл только для функций ode15s,23s,23t,23tb. Если в однойиз этих функций параметр Jacobian не задан, матрица-якобиан строится численно, чтозамедляет работу.Функции ode15s и ode23t можно также применять для решения дифференциальноалгебраических систем.
Они возникают, если известен хотя бы один первый интегралсистемы дифференциальных уравнений, т.е. соотношение между неизвестными функциями,которое обращается в тождество при подстановке в него любого решения системы. Вфизике такие первые интегралы обычно называются законами сохранения (энергии, импульсаи т.д.).Рассмотрим систему дифференциальных уравнений движения центра масс вгравитационном полеx1 = x 2 ,x 2 =-g.Очевидным следствием этой системы является соотношение x 2 x 2 +g x1 =0. Еслипроинтегрировать его, получится 0,5 x22 +g x1 =C.
Это закон сохранения энергии (первоеслагаемое – кинетическая энергия материальной точки с массой единица, второе слагаемое –ее потенциальная энергия). Константа выражается через начальные условия: C=0,5 x22 (0)+g x1 (0). Заменив одно из дифференциальных уравнений первым интегралом, получимдифференциально-алгебраическую системуx1 = x 2 ,0,5 x22 +g x1 -C=0.Тогда функция, вычисляющая правые части, примет следующий вид:function dxdt=funa(t,x);dxdt=[x(2); (1/2)*(x(2)2+ g* x(1) -C];Вторая компонента выходного вектора этой функции в процессе интегрированиядолжна оставаться равной нулю.
Чтобы отразить это обстоятельство, введем матрицумасс M=[1 0; 0 0], и укажем, что эта матрица вырожденная. Оба условия зададим командойoptions=odeset('Mass',M,’MassSingular’,Yes);(t,x)=ode23t(@funa,tspan,x0,options).ПРИЛОЖЕНИЕ 2Решение краевых задач в MATLAB методом стрельбыРассмотрим метод стрельбы на примере решения уравнения второго порядка,предварительно представив его в виде системы уравнений первого порядка.Метод стрельбы базируется на том, что имеются удобные способы численногорешения задачи Коши, т. е. задачи следующего видаy'' = f(x,y,y’), y(0) = y0, y’(0) = tg,(1)где y0 — ордината точки (0,y0), из которой выходит интегральная кривая (рис.2); — уголнаклона интегральной кривой к оси x при выходе из точки (0,y0)._____________________________________________________________________________ 6Деменков Н.П. Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________Рис.2При фиксированном y0 решение задачи (1) имеет вид y=y(x.).
При x=xk (на рис.2 xk=1)решение зависит только от :y( x, ) | x1 =y(1,).Используя указанное замечание о решении задачи Коши, о задачу (1) переформулируемследующим образом: найти такой угол =*, при котором интегральная кривая, выходящаяиз точки (0,Y0) под углом к оси абсцисс, попадет в точку (1,Y1):y(1,)=Y1.(2)Решение задачи (1) при этом =* совпадает с искомым решением задачи (рис.3), т.е.сводится к решению уравнения (2) .Рис.3Уравнение (2) — это уравнение видаF()=0,(3)где F()=y(1,)-Y1.Оно отличается от привычных уравнений лишь тем, что функция F() задана неаналитическим выражением, а с помощью алгоритма численного решения задачи (2).Для решения уравнения (3) можно использовать любой метод, пригодный для уточнениякорней нелинейного уравнения, например, метод деления отрезка пополам, метод Ньютона(касательных) и др.Если имеется достаточно хорошее начальное приближение, метод Ньютонапредпочтительнее из-за высокой стоимости вычисления одного значения функции F(α)(нужно решить задачу Коши (2) с данным ).Метод стрельбы хорошо работает в том случае, если решение y(x,) «не слишкомсильно» зависит от α.