Часть трудности при решении некоторых систем ОДУ определяет подходящее время, чтобы остановить решение. Итоговое время в интервале интегрирования может быть задано определенным событием а не номером. Примером является яблоко, падающее от дерева. Решатель ОДУ должен остановиться, если яблоко поражает землю, но вы не можете знать, когда то событие имело бы место заранее. Точно так же некоторые проблемы включают события, которые не отключают решение. Примером является луна, вращающаяся вокруг планеты. В этом случае вы не можете хотеть останавливать интегрирование рано, но вы все еще хотите обнаружить каждый раз, когда луна завершает один период вокруг большего тела.
Используйте функции события, чтобы обнаружить, когда определенные события будут иметь место во время решения ОДУ. Функции события берут выражение, которое вы задаете и обнаруживаете событие, когда то выражение равно нулю. Они могут также сигнализировать, чтобы решатель ОДУ остановил интегрирование, когда они обнаруживают событие.
Используйте опцию 'Events' функции odeset, чтобы задать функцию события. Функция события должна иметь общую форму
[value,isterminal,direction] = myEventsFcn(t,y)
В случае ode15i функция события должна также принять третий входной параметр для yp.
value выходных аргументов, isterminal и direction являются векторами, i которых th элемент соответствует i th событие:
value(i) является математическим выражением, описывающим i th событие. Событие имеет место, когда value(i) равен нулю.
isterminal(i) = 1, если интегрирование должно остановиться, когда i th событие происходит. В противном случае это - 0.
direction(i) = 0, если все нули должны быть расположены (значение по умолчанию). Значение +1 определяет местоположение только нулей, где функция события увеличивается, и -1 определяет местоположение только нулей, где функция события уменьшается. Задайте direction = [], чтобы использовать значение по умолчанию 0 для всех событий.
Снова, рассмотрите случай яблока, падающего от дерева. ОДУ, которое представляет падающее тело,
с начальными условиями и . Можно использовать функцию события, чтобы определить когда , который является, когда яблоко поражает землю. Для этой проблемы функция события, которая обнаруживает, когда яблоко поражает землю,
function [position,isterminal,direction] = appleEventsFcn(t,y) position = y(1); % The value that we want to be zero isterminal = 1; % Halt integration direction = 0; % The zero can be approached from either direction
Если вы задаете функцию событий, то вызовите решатель ОДУ с тремя дополнительными выходными аргументами, как
[t,y,te,ye,ie] = odeXY(odefun,tspan,y0,options)
Три дополнительных выходных параметра, возвращенные решателем, соответствуют обнаруженным событиям:
te является вектор-столбцом времен, в которые события имели место.
ye содержит значение решения в каждое из времен события в te.
ie содержит индексы в вектор, возвращенный функцией события. Значения указывают, какое событие решатель обнаружил.
Также можно вызвать решатель с одним выводом, как
sol = odeXY(odefun,tspan,y0,options)
В этом случае информация о событии хранится в структуре как sol.te, sol.ye и sol.ie.
Находящий корень механизм, используемый решателем ОДУ в сочетании с функцией события, имеет эти ограничения:
Если терминальное событие имеет место во время первого шага интегрирования, то решатель указывает событие как нетерминальное и продолжает объединяться.
Если больше чем одно терминальное событие имеет место во время первого шага, то только первые регистры события и решатель продолжают объединяться.
Нули определяются пересечениями знака между шагами. Поэтому нули функций с четным числом пересечений между шагами могут быть пропущены.
Если прошедшие события шагов решателя, попытайтесь уменьшать RelTol и AbsTol, чтобы улучшить точность. Также установите MaxStep помещать верхнюю границу в размер шага. Корректировка tspan не изменяет шаги, сделанные решателем.
Этот пример показывает, как записать простую функцию события для использования с решателем ОДУ. Модели ballode файла в качестве примера движение прыгающего мяча. Функция событий останавливает интегрирование каждый раз возвраты шара, и интегрирование затем перезапускает с новыми начальными условиями. Когда шар возвращается, остановки интегрирования и перезапуски несколько раз.
Уравнения для прыгающего мяча
Возврат шара происходит, когда высота шара
равна нулю после уменьшения. Функция событий, которая кодирует это поведение,
function [value,isterminal,direction] = bounceEvents(t,y) value = y(1); % Detect height = 0 isterminal = 1; % Stop the integration direction = -1; % Negative direction only
Введите ballode, чтобы запустить пример и проиллюстрировать использование функции событий, чтобы моделировать возврат шара.
ballode
Этот пример показывает, как использовать направленные компоненты функции события. orbitode файла в качестве примера моделирует ограниченную задачу трех тел, где одно тело вращается вокруг двух намного больших тел. Функция событий определяет точки в орбите, где орбитальное тело является самым близким и самым дальним далеко. Поскольку значение функции событий является тем же самым в самых близких и самых дальних точках орбиты, направление нулевого пересечения - то, что отличает их.
Уравнения для ограниченной задачи трех тел
где
Первые два компонента решения являются координатами тела бесконечно малой массы, таким образом строение один против другого дает орбиту тела.
Функция событий, вложенная в orbitode.m, ищет два события. Одно событие определяет местоположение точки максимального расстояния от отправной точки, и другой определяет местоположение точки, где космический корабль возвращается к отправному вопросу. События расположены точно, даже при том, что размеры шага, используемые интегратором, не определяются местоположениями событий. В этом примере способность задать направление нулевого пересечения очень важна. И точка возвращается к отправному вопросу и точке максимального расстояния от отправной точки, имеют те же значения события, и направление пересечения используется, чтобы отличить их. Функция событий, которая кодирует это поведение,
function [value,isterminal,direction] = orbitEvents(t,y) % dDSQdt is the derivative of the equation for current distance. Local % minimum/maximum occurs when this value is zero. dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4)); value = [dDSQdt; dDSQdt]; isterminal = [1; 0]; % stop at local minimum direction = [1; -1]; % [local minimum, local maximum] end
Введите orbitode, чтобы запустить пример.
orbitode
This is an example of event location where the ability to specify the direction of the zero crossing is critical. Both the point of return to the initial point and the point of maximum distance have the same event function value, and the direction of the crossing is used to distinguish them. Calling ODE45 with event functions active... Note that the step sizes used by the integrator are NOT determined by the location of the events, and the events are still located accurately.
