Неотрицательное решение для ОДУ

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

Если определенные компоненты решения должны быть неотрицательными, то используйте odeset установить NonNegative опция для индексов этих компонентов. Эта опция не доступна для ode23sode15i, или для неявных решателей (ode15sode23tode23tb) примененный проблемы с большой матрицей. В частности, вы не можете наложить ограничения неотрицательности на проблему ДАУ, которая обязательно имеет сингулярную большую матрицу.

Пример: Функция абсолютного значения

Рассмотрите задачу с начальными значениями

$y^{'} = - | y |,$

решенный на интервале $[0, 40]$ с начальным условием $y (0) = 1$ . Решение этого ОДУ затухает, чтобы обнулить. Если решатель производит отрицательное значение решения, то он начинает отслеживать решение ОДУ через это значение, и расчет в конечном счете перестал работать, когда расчетное решение отличается к $- \infty$ . Используя NonNegative опция предотвращает этот отказ интегрирования.

Сравните аналитическое решение $y (t) = e^{- t}$ к решению ОДУ с помощью ode45 без дополнительных опций, и одной с NonNegative опция установлена.

ode = @(t,y) -abs(y);

% Standard solution with |ode45|
options1 = odeset('Refine',1);
[t0,y0] = ode45(ode,[0 40],1,options1);

% Solution with nonnegative constraint
options2 = odeset(options1,'NonNegative',1);
[t1,y1] = ode45(ode,[0 40],1,options2);

% Analytic solution
t = linspace(0,40,1000);
y = exp(-t);

Постройте эти три решения для сравнения. Внушительная неотрицательность крайне важна, чтобы помешать решению повернуть прочь к $- \infty$ .

plot(t,y,'b-',t0,y0,'ro',t1,y1,'k*');
legend('Exact solution','No constraints','Nonnegativity', ...
       'Location','SouthWest')

Пример: Проблема колена

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

$ϵ y^{'} = (1 - x) y - y^{2},$

решенный на интервале $[0, 2]$ с начальным условием $y (0) = 1$ . Параметр $ϵ$ обычно берется, чтобы удовлетворить $0 < ϵ ≪ 1$ , и эта проблема использование $ϵ = 1 \times 1 0^{- 6}$ . Решение этого ОДУ подходы $y = 1 - x$ для $x < 1$ и $y = 0$ для $x > 1$ . Однако вычисление числового решения с допусками по умолчанию показывает, что решение следует $y = 1 - x$ изоклина для целого интервала интегрирования. Внушительные ограничения неотрицательности приводят к правильному решению.

Решите задачу колена с и без ограничений неотрицательности.

epsilon = 1e-6;
y0 = 1;
xspan = [0 2];
odefcn = @(x,y,epsilon) ((1-x)*y - y^2)/epsilon;

% Solve without imposing constraints
[x1,y1] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0);

% Impose a nonnegativity constraint
options = odeset('NonNegative',1);
[x2,y2] = ode15s(@(x,y) odefcn(x,y,epsilon), xspan, y0, options);

Постройте решения для сравнения.

plot(x1,y1,'ro',x2,y2,'b*')
axis([0,2,-1,1])
title('The "knee problem"')
legend('No constraints','Non-negativity')
xlabel('x')
ylabel('y')

Ссылки

[1] Шемпин, L.F., С. Томпсон, Дж.А. Кирженка и Г.Д. Бирн, "Неотрицательные решения ОДУ", Прикладная математика и Издание 170, 2005 Расчета, стр 556-569.

Смотрите также

odeset

Документация