Простые ОДУ, которые имеют компонент единого решения, могут быть заданы как анонимная функция в вызове решателя. Анонимная функция должна принять два входных параметров (t,y), даже если одни из входных параметров не используются в функции.

Решите ОДУ

Задайте временной интервал [0 5] и начальное условие y0 = 0.

tspan = [0 5];
y0 = 0;
[t,y] = ode89(@(t,y) 2*t, tspan, y0);

Постройте решение.

plot(t,y,'-o')

Уравнение Ван дер Поля является ОДУ второго порядка

где $$\mu >0$$ скалярный параметр. Перепишите это уравнение как систему ОДУ первого порядка путем создания замены $$y_1^{\prime }=y_2$$. Получившаяся система ОДУ первого порядка

Файл функции vdp1.m представляет использование уравнения Ван дер Поля $$\mu =1$$. Переменные $$y_1$$ и $$y_2$$ записи y(1) и y(2) из двухэлементного векторного dydt.

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Решите ОДУ с помощью ode89 функция на временном интервале [0 20] с начальными значениями [2 0]. Получившийся выход является вектор-столбцом моментов времени t и массив решения y. Каждая строка в y соответствует времени, возвращенному в соответствующей строке t. Первый столбец y соответствует $$y_1$$, и второй столбец соответствует $$y_2$$.

[t,y] = ode89(@vdp1,[0 20],[2; 0]);

Постройте решения для $$y_1$$ и $$y_2$$ против t.

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE89');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

ode89 работает только с функциями, которые используют два входных параметра, t и y. Однако можно передать дополнительные параметры путем определения их вне функции и передачи их в том, когда вы задаете указатель на функцию.

Решите ОДУ

При перезаписи уравнения, когда уступает система первого порядка

odefcn, локальная функция в конце этого примера, представляет эту систему уравнений как функцию, которая принимает четыре входных параметра: tYA, и B.

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

Решите ОДУ с помощью ode89. Задайте указатель на функцию так, чтобы он передал предопределенные значения для A и B к odefcn.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode89(@(t,y) odefcn(t,y,A,B), tspan, y0);

Постройте график результатов.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

По сравнению с ode45, ode113, ode78, и ode89 решатели лучше в решении задач со строгими ошибочными допусками. Общая ситуация, где эти решатели выделяются, находится в орбитальных проблемах динамики, где кривая решения является гладкой и требует высокой точности на каждом шаге от решателя.

Проблема 2D тела считает две взаимодействующих массы m1 и m2 перемещения по кругу вокруг общей плоскости. В этом примере одна из масс значительно больше, чем другой. С тяжелым телом в начале координат уравнения движения

где

Чтобы решить задачу, сначала преобразуйте в систему четырех ОДУ первого порядка с помощью замен

Замены создают систему первого порядка

twobodyode, локальная функция, включенная в конце этого примера, кодирует систему уравнений для проблемы 2D тела.

function dy = twobodyode(t,y)
% Two-body problem with one mass much larger than the other.
r = sqrt(y(1)^2 + y(3)^2);
dy = [y(2); 
     -y(1)/r^3;
      y(4);
     -y(3)/r^3];
end

Решите систему ОДУ с помощью ode89. Задайте строгие ошибочные допуски 1e-13 для RelTol и 1e-14 для AbsTol.

opts = odeset('Reltol',1e-13,'AbsTol',1e-14,'Stats','on');
tspan = [0 10*pi];
y0 = [2 0 0 0.5];
[t,y] = ode89(@twobodyode, tspan, y0, opts);

243 successful steps
0 failed attempts
5103 function evaluations

plot(t,y)
legend('x','x''','y','y''','Location','SouthEast')
title('Position and Velocity Components')

figure
plot(y(:,1),y(:,3),'-o',0,0,'ro')
axis equal
title('Orbit of Smaller Mass')

По сравнению с ode45, ode89 решатель может получить решение быстрее и с меньшим количеством шагов и вычислений функции.

function dy = twobodyode(t,y)
% Two-body problem with one mass much larger than the other.
r = sqrt(y(1)^2 + y(3)^2);
dy = [y(2); 
    -y(1)/r^3;
    y(4);
    -y(3)/r^3];
end