Решение простой системы дифференциальных уравнений - PullRequest
Новая
       101

Решение простой системы дифференциальных уравнений

0 голосов
/ 18 февраля 2020

Как бы я численно решил следующую простую систему дифференциальных уравнений, используя Октава ?

enter image description here

enter image description here

Примечание:

  • Я использую квалификатор «простой», так как, насколько я понимаю, система первого порядка и не связана .
  • Я попробовал все методы и сценарии онлайн, чтобы попытаться решить эту проблему, включая здесь , здесь и здесь . Во всех вариантах я получаю либо зависшую, не реагирующую Octave, подсказку о «повторяющихся сбоях конвергенции», ошибку с рекомендацией вручную отрегулировать начальный и максимальный размер шага (что я пробовал и делал, но безрезультатно) или что-то, что изначально кажется решением из-за отсутствия ошибок, но на графике решения показан пустой график
  • Где Octave предоставил эквивалентные подпрограммы Matlab, я пробовал различные подпрограммы ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, ode15i и, конечно же, октавы владеют командой lsode, и все они дают те же ошибки, что и описанные выше.

1 Ответ

4 голосов
/ 18 февраля 2020

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

% z = [x,y]
f = @(t,z) [ z(1).^2+t; z(1).*z(2)-2 ];
z0 = [ 2; 1];

[ T, Z ] = ode45(f, [0, 10], z0);

plot(T,Z); legend(["x";"y"]);

Интегратор не работает, о чем сообщается с предупреждением

: решение не выполнено. Итеративная интеграция l oop вышла за время t = 0.494898 до того, как была достигнута конечная точка на tend = 10.000000. Это может произойти, если размер шага становится слишком маленьким. Попробуйте уменьшить значение 'InitialStep' и / или 'MaxStep' с помощью команды 'odeset'.

Повторение интегрирования незадолго до критического времени 

opt = odeset('MaxStep',0.01);
[ T, Z ] = ode45(f, [0, 0.49], z0, opt);
clf; plot(T,Z); legend(["x";"y"]);

приводит к графику

plot of the solution

, где можно видеть, что квадратичный член c в первом уравнении приводит к убеганию. По какой-то причине решатель распознает только уменьшающийся размер шага, но не убегающие значения решения.

Действительно, первое - это уравнение Риккати, которое, как известно, имеет полюса в конечное время. Используя типичную параметризацию, x(t)=-u'(t)/u(t) имеет по правилу произведение / частное производную 

x' = -u''(t)/u(t) - u'(t)* (-u'(t)/u(t)^2) = -u''(t)/u(t) + x(t)^2

, которая затем приводит к ODE для u

u''(t)+t*u(t)=0, u(0)=-1, u'(0)=x(0)=2,

, которое является уравнением Эйри с Колеблющаяся ветвь для t>0. Первый root из u является полюсом для x, нет возможности расширить решение за эту точку. 

g=@(t,u) [u(2); -t.*u(1)]
u0 = [ 1; -2];

function [val,term, dir] = event(t,u)
   val = u(1);
   term = 0;
   dir = 0;
end
opt = odeset('MaxStep',0.1, 'Events', @(t,u) event(t,u));
[T,U,Te,Ue,Ie] = ode45(g,[0,4],u0,opt);
disp(Te)
clf; plot(T,U); legend(["u";"u'"]);

, в котором нули u перечислены как 0.4949319379979706, 2.886092605590324 , еще раз подтверждая причину предупреждения, и выдает сюжет

enter image description here

...