Случайное событие Matlab/Octave ode45

У меня есть некоторые проблемы с пониманием того, как реализовать события в октаве/матлабе при решении дифференциальных уравнений.

Рассмотрим, например, этот простой код для решения дифференциального уравнения y' = -y:

function dy = odefun(t,y)
    dy = -y;
endfunction

options = odeset('RelTol',1e-4,'AbsTol',[1e-4]);
[T,Y] = ode45(@odefun,[0 12],[1],options);

Теперь я хотел бы представить случайное событие. Например, на каком-то фиксированном временном шаге я хотел бы случайным образом изменить значение y, а затем продолжить эволюцию в соответствии с дифференциальным уравнением. Как я могу это сделать?


person Leonardo    schedule 01.02.2014    source источник
comment
Кажется, это дубликат этого вопроса (ваш вопрос намного лучше), но он не может быть помечен как таковой, потому что ни один из ответов не был принят. Поэтому я собираюсь опубликовать свой ответ здесь и настроить его на то, что вы просите.   -  person horchler    schedule 01.02.2014


Ответы (1)


Способ сделать это состоит в том, чтобы интегрировать систему по частям и объединять результирующие выходные данные каждого прогона:

% t = 0 to t = 6
tspan = [0 6];
y0 = 1;
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[T,Y] = ode45(@odefun,tspan,y0,options);

% t = 6 to t = 12
tspan = [6 12];
y0 = Y(end,:)+randn; % Pertubation
[t,y] = ode45(@odefun,tspan,y0,options);
T = [T;t(2:end)];  % Remove first value as it will be same as last of previous run
Y = [Y;y(2:end,:)];

Редактор Matlab может пожаловаться на то, что массивы T и Y не были предварительно выделены и/или не росли, но в этом случае это нормально, так как они растут большими кусками всего несколько раз.

Чего вы не хотите пытаться сделать (но многие пытаются), так это добавить свое возмущение внутрь функции интегрирования. Во-первых, если вы хотите, чтобы это произошло точно в определенное время или при соблюдении определенных условий, это нельзя сделать из функции ОДУ. Во-вторых, вставка больших разрывов в ОДУ может привести к снижению точности и увеличению времени вычислений (особенно с ode45 - ode15s может быть лучшим вариантом или, по крайней мере, убедитесь, что ваши абсолютные и относительные допуски подходят). Вы бы фактически создали очень жесткую систему.

Если вы хотите, чтобы ваше возмущение возникало при выполнении определенных условий (а не в определенное время), вам нужно будет использовать функцию события (см. этот вопрос, например). Вам нужно будет использовать функцию событий, чтобы завершить интегрирование, когда условия будут выполнены, а затем применить свое возмущение, как я показал выше.

person horchler    schedule 01.02.2014