Простые ОДУ, имеющие один компонент решения, могут быть указаны как анонимная функция в вызове решателя. Анонимная функция должна принимать два входа (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 решатель - хороший первый выбор для большинства жестких задач. Однако другие жесткие решатели могут быть более эффективными для определенных типов проблем. В этом примере решается жесткое тестовое уравнение с использованием всех четырех жестких решателей ОДУ.

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

Уравнение становится все более жестким по мере $$$$увеличения величины λ. Используйте $$\lambda \mathrm{=}1{\mathrm{}}^{\mathrm{\times }}$$ 109 и исходное $$\mathrm{условие}y\mathrm{}$$(0) = 1 в течение интервала времени[0 0.5]. С этими значениями проблема достаточно жесткая, что ode45 и ode23 бороться за интеграцию уравнения. Кроме того, используйте odeset для передачи постоянного $$\frac{}{}\mathrm{J=\partial f\partial y=-\lambda }$$ Якобиана и включения отображения статистики решателя.

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.338321 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.258983 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.406110 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.355101 seconds.

title('ode23tb')

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

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

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

Решите уравнение van der Pol с $$\lambda \mathrm{=}\mathrm{}\mathrm{}\mathrm{}$$1000 с помощью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 для генерации 2500 точек в интервале [0 3000]. Вычислите первый компонент решения в этих точках с помощью deval.

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

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

plot(x,y)

Расширить решение до $${}_{\mathrm{tf}}\mathrm{=}\mathrm{}\mathrm{}\mathrm{}$$4000 с помощью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')

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

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

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

Система уравнений может быть переписана как система DAE, используя закон сохранения для определения состояния. При этом проблема преобразуется в систему дисковой полки.

Дифференциальный индекс этой системы равен 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');