Чтобы найти численное решение в matlab, нужно определить
VDP_ODE = @(t,u) [ u(2), 2*a*u(2)*(1-u(1)^2) ]; % u = [ y, z ]
(или определить его с помощью полной декларации function) и вызвать его с помощью числового решателя
a=0.025; h=0.1;
t = 0:h:10;
y0 = 0; z0 = 0.5;
u0 = [y0, z0 ];
u = ode45(@VDP_ODE, t, u0)
figure(1); plot(t,u(:,1));
figure(2); plot(u(:,1),u(:,2));
Если вы хотите создать свой собственный решатель, вы должны сначала понять метод.Кажется, что это метод 4-го порядка, подобный, но не совсем равный линейному многоступенчатому методу, подобно методам Бимана-Шофилда, введенным для вычисления молекулярной динамики.
Входными данными для каждого шага являются значения y(n), y(n+1), выходные данные уменьшаются до y(n+2) без других значений, перенесенных на следующий шаг.Внутри шага метода дополнительные значения вычисляются для y(n+3/2) и производных z во все времена выборки t(n), t(n+1), t(n+3/2), t(n+2).Цель состоит в том, чтобы получить результаты с плавающей запятой, поэтому нет смысла определять уравнения символически, тем самым адаптировать систему к интерфейсу fsolve.
Система, которая решается в коде Maple, может быть упрощенанаблюдая, что есть 6 уравнений, таким образом, 6 неизвестных параметров.Метод основан на полиноме 5-й степени с 6 коэффициентами в качестве неизвестных, который интерполирует значения, первую и вторую производные решения.Это дает 2 уравнения для заданных значений y0,y1 и 4 уравнения для точного удовлетворения дифференциального уравнения в точках интерполяции.
Y = [ y0, step0(y0,z0) ]
for n=1:N-2
Y(n+2,:)=stepN(Y(n,:), Y(n+1,:), h)
end
и для шага метода сделайте что-то вроде
%%%% ======== General step ========
%% fit a degree 5 polynomial to values, 1st and 2nd derivatives
%% at t0, t1, t1h=t1+0.5*h, t2
%% p(s) = sum c_k s^k/k!, c0=y1, y(t1+s) ~ p(s)
%% eqn y0 = p(-h)
%% def z0 = p'(-h)
%% eqn f(y0,z0) = p''(-h)
%% def z1 = p'(0)
%% eqn f(y1,z1) = p''(0)
%% eqn y2 = p(h), etc
%%
%% state vector for non-lin solver is u = [y2, c1, c2, ...c5 ]
function y2 = stepN(y0,y1,h)
zz = (y1-y0)/h;
predict = [ y1+h*zz, zz, f(y1,zz), 0, 0, 0];
options = optimset("TolX", 1e-1*h^6, "TolFun", 1e-2*h^6);
u = fsolve(@(u) systemN(u,[y0,y1], h), predict, options);
y2 = u(1);
end
, который устанавливает некоторые значения предиктора (только O (h)), устанавливает допуски для решателя, так что ошибка решателя, как мы надеемся, меньше ошибки дискретизации, равной O(h^6), вызывает конец решателя, извлекает искомыйзначение из массива решения.
Чтобы система могла решить функцию, сначала извлекает константы и переменные в именованные переменные для лучшей читаемости, вычисляет значения функции в данной точке, а затем возвращает массив дефектов в дискретизированных уравнениях.
function eqn = systemN(u,init,h)
[ y0, y1 ] = num2cell(init){:};
[ y2, c1, c2, c3, c4, c5 ] = num2cell(u){:};
p = @(h) y1 + h*(c1+h/2*(c2+h/3*(c3+h/4*(c4+h/5*c5))));
dp = @(h) c1+h*(c2+h/2*(c3+h/3*(c4+h/4*c5)));
d2p = @(h) c2+h*(c3+h/2*(c4+h/3*c5));
z0 = dp(-h); z1 = c1;
y1h = p(0.5*h); z1h = dp(0.5*h);
z2 = dp(h);
eqn = [ y2-p(h),
y0-p(-h),
f(y0,z0)-d2p(-h),
f(y1,z1)-c2,
f(y1h,z1h)-d2p(0.5*h),
f(y2,z2)-d2p(h) ];
end