Численные методы ОДУ. Чувствительность от шага

[robot80]

Здравствуйте, все!
Тестируя собственноручно написанный явный метод Адамса для ОДУ первого порядка, встретился с весьма неприяной ситуацией (и для всех стандартных основных методов, кстати тоже) зависимости точности решения от (наверное, погрешности определения) величины шага. Код метода в scilab вот:

function [tr, xr] = Adams4 (t0,x0,h,te)
tr = [t0: h: te]; n = length(tr);xr =zeros (n); 
[t1, x1] = RK4 (t0, x0, h); [t2, x2] = RK4 (t1, x1, h);
[t3, x3] = RK4 (t2, x2, h);
xr(1) = x0; xr(2) = x1;xr(3) = x2; xr(4) = x3;
f0 = func(tr(1), xr(1)); f1 = func(tr(2), xr(2));
f2 = func(tr(3), xr(3)); f3 = func(tr(4), xr(4));
for k = 5:1:length(tr)
dx=55.0.*f3-59.0.*f2+37.0.*f1-9.0.*f0;
xr(k) = xr(k-1)+ h.*dx./24.0; 
f0 = f1; f1 = f2; f2 = f3; f3 = func(tr(k), xr(k));
end
endfunction

function [t,x] = RK4 (t,x,h, t0,te)
g1 = h.*func(t,x);
t2 = t + 0.5.*h; x2 = x + 0.5.*g1;
g2 = h.*func(t2, x2);
t3 = t + 0.5.*h; x3 = x + 0.5*g2;
g3 = h.*func(t3, x3);
t4 = t + h; x4 = x + g3;
g4 = h.*func(t4, x4);
dx=(g1+2.0.*g2+2.0.*g3+ g4)./6.0;
x = x + dx; t = t+h;
endfunction

function res = func(t,x)
res = -x.*tan(t) + 1.0./cos(t);
endfunction

Так вот, делая два следующих вызова этой функции:

[te, xe2] = Adams4 (%pi, 1.0, 3.0.*%pi./99.0, 4.0.*%pi);

и

[te, xe2] = Adams4 (%pi, 1.0, 3.0.*%pi./100.0, 4.0.*%pi);

, получаем в первом случае правильный ответ, а во втором неправильный совершенно.
Любопытно, что встроенный в scilab решатель в обоих случаях решает правильно.
И вот по по этому поводу вопрос: можно ли с этой сильной неустойчивостью чего-нибудь сделать, как-нибудь её победить?
Как гипотеза: использовать переменный шаг при этом подбирать его так, чтобы он согласовался с основным шагом.

[Pphantom]

robot80 в сообщении #1331729 писал(а):

получаем в первом случае правильный ответ, а во втором неправильный совершенно

А что Вы считаете в данном случае "правильным ответом"? На интервале интегрирования $t \in [\pi, 4\pi]$ правая часть очевидно разрывна.

[Евгений Машеров]

Изменение величины шага на результат решения повлияет, но не изменение на 1%.
Могу только предположить, что функция разрывная, и при перемене шага то ли попадаем на точки разрыва, то ли не попадаем. И, собственно, проблема лишь в том, что метод, предполагающий непрерывность функции, не обязательно работает на разрывных (и необязательно не работает, "как повезёт").

[robot80]

Pphantom, Правильный ответ в данном случае тот, который максимально совпадает с аналитическим решением.
Евгений Машеров, да, похоже, действительно, узел сетки попадал в разрыв функции. Удивительно, конечно, такое точное совпадение двух чисел. Я как-то сразу об этом не подумал. А внутренний сайлабовский решатель, видимо, как-то их обходит.
Спасибо за совет! Тему можно закрыть.

[Pphantom]

robot80 в сообщении #1331772 писал(а):

Pphantom, Правильный ответ в данном случае тот, который максимально совпадает с аналитическим решением.

В таком случае есть проблема. Продолжение решения задачи Коши за точку разрыва само по себе некорректно, для соответствующего деяния само понятие решения нужно некоторым образом обобщать. Пока такого обобщения нет, любой численный результат для $t>3 \pi/2$ в любом случае неверен (даже если внешне соответствует интуитивным представлениям о "правильности"). Непосредственно при счете, конечно, проблема возникает по указанной Евгений Машеров причине, но еще вопрос, что лучше - получить очевидно нелепый ответ или случайно более или менее удачно "проскочить" точку разрыва. Во втором случае ответ может показаться правдоподобным (и при этом останется неверным даже в рамках тех самых интуитивных представлений).

[DeBill]

Во втором примере, пятидесятоя точка точно вляпывается в особенность. Но, поскольку косинус от пи плюс стопятьдесят сотых пи немножко не равен нулю, программа не вылетает (а просто выдает чушь)...
Забавно, как бы эту задачу (найти значение в точку 4пи) решали:
а) первокурсник: запускает Вольфрам; в ответ $x=\sin t -\cos t$ подставляет тэ, получает -1
б) ленивый прикладник: запускает пакет, получает -1, но иногда (по пятницам, 13) 10000;
б1) неленивый: реализует ручками какойнить численный метод, с тем же результатом - но по всем пятницам
Более опытный прикладник объясняет им, что - при взбрыке программы - они нарушили Священный Ритуал: при выборе шага, надо, чтоб коэффициент при пи содержал в знаменателе простое число, большее 17. А при запуске программы, надо в левой руке держать банку с пивом...
в) второкурсник: решает линейное неоднородное ур-е, далее как а)
г) старшекурсник с дифуров, качественная теория: находит особую точку $x=-1, t= \frac{3\pi}{2}$ , ее орбитальный тип - "устойчивый узел" , делает вывод, что: все траектории вляпываются в эту точку, дальше нет смысла ее продолжать.
д) старшекурсник с дифуров, гладкая теория: замечает, что узел - дикритический, так что траектории гладко продолжаются. Делает вывод: решение в точке 4пи можно найти - численно, если повезет.
е) препод: выходит в комплексную область, исследует все особые точки римановой поверхности "решение", убеждается что все они - однозначного характера, аналитически продолжает решение вокруг всех особых точек, и делает вывод: решение в точке 4пи существует...