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

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

решенный на интервале $$[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. Уравнение

решенный на интервале $$[0,2]$$ с начальным условием $$y(0)=1$$. Параметр $$ϵ$$ обычно берется, чтобы удовлетворить $$0<ϵ\ll 1$$, и эта проблема использование $$ϵ=1\times 10^{-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.