Matlab: найти корень значений y, заданных дифференциальным уравнением

Question

Я решил начальное дифференциальное уравнение и построил графики значений Y, заданных ODE45. По сюжету я смутно могу сказать, где должен быть корень, но в данной задаче мне нужно найти его с большой точностью.

Моим первым предположением было настроить полином на мои значения X и Y, а затем решить уравнение полинома. Но я использовал polyfit и имел 69 значений, которые давали мне полином 68-градусного, который я не мог решить. Итак, кто-нибудь знает, как я могу найти "корень" для набора значений Y, не зная фактического уравнения? В задании написано, что следует использовать интерполяцию!

Заранее спасибо!

Answer 1

Учитывая вектор значений Y (упорядоченный в том смысле, что соответствующие значения X неуклонно растут), вы можете легко найти, возле каких значений X расположены корни. Корни либо там, где значение Y равно нулю, либо между двумя последовательными значениями Y, которые меняют знак. Идея проиллюстрирована в следующем фрагменте кода:

X = -1:0.1:1;
Y = X.*X - 0.4;

root_exact_pos  = find(Y==0);
root_approx_pos = find(diff(sign(Y))~=0);

Корни в X значения, либо в X(root_exact_pos(k))или между X(root_approx_pos(k)) а также X(root_approx_pos(k)+1), k переходя от 1 к числу элементов соответствующего массива корневых позиций.

С этого момента вы можете применять любую интерполяцию, которую хотите, чтобы найти лучшее приближение корня (я бы пошел с линейным, между 2 точками).

Answer 2

Вы сказали, что вам нужно найти рут "с большой точностью". Ответ @CST-Link выше не в том, как это сделать, если вы решаете дифференциальное уравнение численно с помощью ODE45, На самом деле это плохая идея. Это требует вывода точек вашего решения с высоким разрешением и вносит ошибку, которая может быть особенно плохой, если интегрируемые уравнения являются консервативными (т. Е. Вы будете эффективно прибавлять или вычитать энергию из решения, в котором полная энергия всегда должна оставаться постоянной). Вам необходимо использовать функцию обнаружения событий (пересечения нуля) решателей ODE Matlab. Как правило, вы сможете найти ваш рут с точностью, близкой к машинному эпсилону, eps, Принимая во внимание, что метод, данный @CST-Link, даст только точность порядка шага ODE45, который может быть очень большим (для улучшения этого можно использовать интерполяцию, но вы не собираетесь приближаться к eps если вы не используете крошечные размеры шага).

Глядя на помощь Matlab для ODE45функциональность событий может показаться запутанной, поэтому я попытаюсь привести более простой пример с использованием кода, основанного на ballode (help ballode а также edit ballode для получения дополнительной информации), пример, который Matlab включает в себя для демонстрации событий.

function eventsdemo

options = odeset('Events',@efun); % specify name of events function
[t,y,te,ye,ie] = ode45(@f,[0 10],[0;20],options); % integrate
figure
plot(t,y(:,1),'b',te,ye(:,1),'r.')

function dydt = f(t,y)
dydt = [y(2);-9.8]; % differential equation for ballistic motion

function [value,isterminal,direction] = efun(t,y)
value = y(1);
isterminal = 1;
direction = -1;

Примером является просто простое баллистическое движение шара в вертикальном направлении: функция f, Цель состоит в том, чтобы точно определить скорость мяча во время его пересечения с земной плоскостью, y(1) == 0 линия в отрицательном направлении. Если вы обнаружите контакт слишком рано, скорость шара будет меньше реальности, и энергия будет потеряна; слишком поздно и энергия будет введена. value вывод efun определяет уравнение, которое должно равняться нулю для пересечения. Это может быть так просто или так сложно, как необходимо, может зависеть от всех переменных состояния и времени, и вы можете обнаружить несколько пересечений, указав вектор. В вашем случае это звучит так, как будто вы заинтересованы в корнях вдоль оси x, поэтому я полагаю, что вы можете иметь идентичное определение value, Если вам нужен только первый корень или есть только один, то isterminal остановит интеграцию, когда она будет найдена. Наконец, если вы не знаете наклона / градиента вашей функции в корне, вы можете установить direction в ноль. Попробуйте эту функцию события с приведенным выше кодом (ie может использоваться, чтобы сказать, какое из трех событий вызвано в te а также ye выходы):

function [value,isterminal,direction] = efun(t,y)
value = [y(2) y(1)-10 y(1)];
isterminal = [0 0 1];
direction = [-1 0 -1];

Некоторые незначительные предостережения. Вы должны установить конечное время интеграции достаточно большим. Другими словами, событие должно происходить между t0 а также tf (увидеть ballode для схемы итерации вызовов ODE45 если вы не знаете, где со временем будут происходить ваши события). Если вы установите isterminal на ноль, выход t вектор и y Матрица будет обрезана до конца во время терминального события. Наконец, если вы укажете фиксированный выходной размер, например, TSPAN = t0:dt:tfпоследняя точка вывода, вероятно, будет иметь меньший размер шага.