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

Решите ОДУ

Используйте временной интервал [0,2] и начального условия y0 = 1.

tspan = [0 2];
y0 = 1;
[t,y] = ode15s(@(t,y) -10*t, tspan, y0);

Постройте график решения.

plot(t,y,'-o')

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

Система уравнений:

Начальные условия и. Функциональный vdp1000 поставляется с MATLAB® и кодирует уравнения.

function dydt = vdp1000(t,y)
%VDP1000  Evaluate the van der Pol ODEs for mu = 1000.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB.

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

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

Решение этой системы с помощью ode45 с допусками относительной и абсолютной погрешности по умолчанию (1e-3 и 1e-6, соответственно) является чрезвычайно медленным, требуя, чтобы несколько минут решили и построили график решения. ode45 требует, чтобы миллионы временных шагов завершили интегрирование, из-за областей жесткости, где это изо всех сил пытается соответствовать допускам.

Это - график решения, полученного ode45, который занимает много времени, чтобы вычислить. Заметьте огромное количество временных шагов, требуемых проходить через области жесткости.

[t,y] = ode15s(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

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

Решите ОДУ

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

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

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

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

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

Решатель ode15s является хорошим предпочтительным вариантом для самых жестких проблем. Однако другие жесткие решатели могут быть более эффективными для определенных типов проблем. Этот пример решает жесткое тестовое уравнение с помощью всех четырех жестких решателей ОДУ.

Рассмотрите тестовое уравнение

Уравнение становится все более и более жестким как значение увеличений. Используйте и начальное условие по временному интервалу [0 0.5]. С этими значениями проблема достаточно жестка, что ode45 и ode23 изо всех сил пытаются интегрировать уравнение. Кроме того, используйте odeset, чтобы передать в постоянном якобиане и включить отображение статистики решателя.

lambda = 1e9;
y0 = 1;
tspan = [0 0.5];
opts = odeset('Jacobian',-lambda,'Stats','on');

Решите уравнение с ode15s, ode23s, ode23t и ode23tb. Сделайте подграфики для сравнения.

subplot(2,2,1)
tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc

104 successful steps
1 failed attempts
212 function evaluations
0 partial derivatives
21 LU decompositions
210 solutions of linear systems
Elapsed time is 1.941491 seconds.

title('ode15s')
subplot(2,2,2)
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc

63 successful steps
0 failed attempts
191 function evaluations
0 partial derivatives
63 LU decompositions
189 solutions of linear systems
Elapsed time is 0.483929 seconds.

title('ode23s')
subplot(2,2,3)
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc

95 successful steps
0 failed attempts
125 function evaluations
0 partial derivatives
28 LU decompositions
123 solutions of linear systems
Elapsed time is 0.632605 seconds.

title('ode23t')
subplot(2,2,4)
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc

71 successful steps
0 failed attempts
167 function evaluations
0 partial derivatives
23 LU decompositions
236 solutions of linear systems
Elapsed time is 0.357163 seconds.

title('ode23tb')

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

Для общих жестких проблем производительность жестких решателей отличается в зависимости от формата проблемы и заданных опций. Обеспечение якобиевской матрицы или шаблона разреженности всегда повышает эффективность решателя для жестких проблем. Но поскольку жесткие решатели используют якобиан по-другому, улучшение может значительно отличаться. В сущности, если система уравнений является очень большой или должна много раз решаться, то стоит исследовать производительность различных решателей, чтобы минимизировать время выполнения.

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

Решите уравнение Ван дер Поля с использованием ode15s. Функциональный vdp1000.m поставляется с MATLAB® и кодирует уравнения. Задайте единственный вывод, чтобы возвратить структуру, содержащую информацию о решении, таком как точки оценки и решатель.

tspan = [0 3000];
y0 = [2 0];
sol = ode15s(@vdp1000,tspan,y0)

sol = struct with fields:
     solver: 'ode15s'
    extdata: [1x1 struct]
          x: [1x592 double]
          y: [2x592 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

Используйте linspace, чтобы сгенерировать 2 500 точек в интервале [0 3000]. Оцените первый компонент решения в этих точках с помощью deval.

x = linspace(0,3000,2500);
y = deval(sol,x,1);

Постройте график решения.

plot(x,y)

Расширьте решение использования odextend и добавьте результат к исходному графику.

tf = 4000;
sol_new = odextend(sol,@vdp1000,tf);
x = linspace(3000,tf,350);
y = deval(sol_new,x,1);
hold on
plot(x,y,'r')

Этот пример повторно формулирует систему ОДУ как система дифференциальных алгебраических уравнений (ДАУ). Задачей Робертсона, найденной в hb1ode.m, является классическая тестовая задача для программ, которые решают жесткие ОДУ. Система уравнений

hb1ode решает эту систему ОДУ к устойчивому состоянию с начальными условиями, и. Но уравнения также удовлетворяют линейный закон сохранения,

С точки зрения решения и начальных условий, закон сохранения

Система уравнений может быть переписана как система ДАУ при помощи закона сохранения, чтобы определить состояние. Это повторно формулирует проблему как систему ДАУ

Дифференциальный индекс этой системы равняется 1, поскольку только единственная производная требуется, чтобы делать это системой ОДУ. Поэтому никакие дальнейшие преобразования не требуются прежде, чем решить систему.

Функциональный robertsdae кодирует эту систему ДАУ. Сохраните robertsdae.m в своей текущей папке, чтобы запустить пример.

function out = robertsdae(t,y)
out = [-0.04*y(1) + 1e4*y(2).*y(3)
   0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2
   y(1) + y(2) + y(3) - 1 ];

Полный пример кода для этой формулировки задачи Робертсона доступен в hb1dae. m.

y0 = [1; 0; 0];
tspan = [0 4*logspace(-6,6)];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]);
[t,y] = ode15s(@robertsdae,tspan,y0,options);

y(:,2) = 1e4*y(:,2);
semilogx(t,y);
ylabel('1e4 * y(:,2)');
title('Robertson DAE problem with a Conservation Law, solved by ODE15S');