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

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

решен на интервале $$[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., S. Thompson, J.A. Kierzenka, and G.D. Byrne, «Non-негативные решения ОДУ», Applied Mathematics and Computation Vol. 170, 2005, pp. 556-569.