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

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

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

Уравнение становится все более жестким как величина $$\lambda$$ увеличивается. Использовать $$\lambda =1\times 10^9$$ и начальное условие $$y(0)=1$$ на протяжении временного интервала [0 0.5]. С этими значениями задача достаточно жесткая, чтобы ode45 и ode23 бороться за интегрирование уравнения. Кроме того, используйте odeset пройти в постоянном якобиане $$J=\frac{\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 наиболее поскольку он обычно оценивает якобиан на каждом шагу.

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

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

Решите уравнение Ван дер Поля с $$\mu =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)

Расширение решения на $$t_f=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')

Этот пример переформулирует систему ОДУ как систему дифференциальных алгебраических уравнений (ДАУ). Задача Робертсона, найденная в 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');