Отчасти трудность в решении некоторых систем ОДУ заключается в определении подходящего времени для остановки решения. Конечное время в интервале интеграции может определяться конкретным событием, а не числом. Пример - яблоко, падающее с дерева. Решатель ODE должен остановиться, как только яблоко ударится о землю, но вы можете не знать, когда это событие произойдет заранее. Аналогичным образом, некоторые проблемы связаны с событиями, которые не завершают решение. Примером может служить луна, вращающаяся вокруг планеты. В этом случае вы, возможно, не хотите остановить интеграцию рано, но вы все еще хотите обнаружить каждый раз, когда Луна завершает один период вокруг большего тела.
Используйте функции событий для обнаружения определенных событий во время решения ОДУ. Функции события принимают указанное выражение и обнаруживают событие, когда это выражение равно нулю. Они также могут сигнализировать решателю ODE остановить интеграцию при обнаружении события.
Используйте 'Events' вариант odeset для указания функции события. Функция события должна иметь общую форму
[value,isterminal,direction] = myEventsFcn(t,y)
В случае ode15i, функция события должна также принять третий входной аргумент для yp.
Выходные аргументы value, isterminal, и direction являются векторами, i-й элемент соответствует i-е событие:
value(i) - математическое выражение, описывающее i-е событие. Событие возникает, когда value(i) равно нулю.
isterminal(i) = 1 если интеграция должна завершиться, когда iПроисходит второе событие. В противном случае это 0.
direction(i) = 0 если должны быть найдены все нули (по умолчанию). Значение +1 находит только нули, где увеличивается функция события, и -1 находит только нули, где функция события уменьшается. Определить direction = [] для использования значения по умолчанию 0 для всех событий.
Рассмотрим случай падения яблока с дерева. ОДУ, который представляет падающее тело
'2,
с начальными условиями = 1 0) = 0. Можно использовать функцию события, чтобы определить, y (t) = 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
Если указана функция событий, вызовите решатель ODE с тремя дополнительными выходными аргументами, как
[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 не изменяет шаги, выполняемые решателем.
В этом примере показано, как записать простую функцию события для использования с решателем ODE. Пример файла 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.
