Exponenta Banner

ode15s

Решение жестких дифференциальных уравнений и ДАУ - метод переменного порядка точности

Свернуть все на странице

Синтаксис

[t,y] = ode15s(odefun,tspan,y0)
[t,y] = ode15s(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode15s(odefun,tspan,y0,options)
sol = ode15s(___)

Описание

пример

[t,y] = ode15s(odefun,tspan,y0), где tspan = [t0 tf], интегрирует систему дифференциальных уравнений y'=f(t,y) от t0 на tf с начальными условиями y0. Каждая строка массива решений y соответствует значению, возвращаемому в вектор-столбец t.

Все MATLAB® Решатели ОДУ могут решить системы уравнений вида y'=f(t,y), или задачи, которые включают большую матрицу, M(t,y)y'=f(t,y). Все решатели используют аналогичные синтаксисы. ode23s решатель может решить задачи только с большой матрицей, если большая матрица постоянна. ode15s и ode23t может решить задачи с большой матрицей, сингулярной, известной как дифференциально-алгебраические уравнения (ДАУ). Задайте большую матрицу используя Mass опция odeset.

пример

[t,y] = ode15s(odefun,tspan,y0,options) также использует настройки интегрирования, заданные как options, который является аргументом, созданным с помощью odeset функция. Для примера используйте AbsTol и RelTol опции для задания абсолютных и относительная погрешность допусков или Mass опция для задания большой матрицы.

[t,y,te,ye,ie] = ode15s(odefun,tspan,y0,options) дополнительно находит, где функции (t, y), называемые функциями события, равны нулю. В выходах te - время события, ye является решением во время события и ie - индекс инициируемого события.

Для каждой функции события задайте, должно ли интегрирование завершаться на нуле и имеет ли значение направление пересечения нуля. Сделайте это, установив 'Events' свойство функции, например myEventFcn или @myEventFcn, и создание соответствующей функции: [value, isterminal, direction] = myEventFcn(t, y). Для получения дополнительной информации смотрите Расположение события ОДУ.

пример

sol = ode15s(___) возвращает структуру, которую можно использовать с deval для оценки решения в любой точке интервала [t0 tf]. Можно использовать любой из комбинаций входных аргументов в предыдущих синтаксисах.

Примеры

свернуть все

Открыть Live Script

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

Решение ОДУ

y′=-10t.

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

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

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

plot(t,y,'-o')

Figure contains an axes. The axes contains an object of type line.

Открыть скрипт

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

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

$$\begin{array}{cl} y_1' &= y_2\\y_2' &= 1000(1-y_1^2)y_2-y_1\end{array}$$

Начальными условиями являются$y_1(0)=2$ и. $y_2(0)=0$Функция 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, вычисление которого занимает много времени. Заметьте огромное количество временных шагов, необходимых для прохождения через области жесткости.

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

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

Открыть скрипт

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

Решение ОДУ

$$y'' = \frac{A}{B} t y.$$

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

$$\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= \frac{A}{B} t y_1.
\end{array}$$

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);

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

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),'-.')

Открыть Live Script

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

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

y′=-λy.

Уравнение становится все более жестким как величина λ увеличивается. Использовать λ=1×109 и начальное условие y(0)=1 на протяжении временного интервала [0 0.5]. С этими значениями задача достаточно жесткая, чтобы ode45 и ode23 бороться за интегрирование уравнения. Кроме того, используйте odeset пройти в постоянном якобиане J=∂f∂y=-λ и включите отображение статистики решателя.

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')

Figure contains 4 axes. Axes 1 with title ode15s contains 2 objects of type line. Axes 2 with title ode23s contains 2 objects of type line. Axes 3 with title ode23t contains 2 objects of type line. Axes 4 with title ode23tb contains 2 objects of type line.

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

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

Открыть Live Script

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

y1′′-μ(1-y12)y1′+y1=0.

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

Figure contains an axes. The axes contains an object of type line.

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

Figure contains an axes. The axes contains 2 objects of type line.

Открыть скрипт

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

$$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 y_2y_3- (3 \times 10^7)y_2^2\\ y'_3 &= (3 \times
10^7)y_2^2.\end{array}$$

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

$$y'_1 + y'_2 + y'_3 = 0.$$

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

$$y_1 + y_2 + y_3 = 1.$$

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

$$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 y_2y_3-(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 - 1.\end{array}$$

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

Функция 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

Решить систему ДАУ можно используя ode15s. Допустимые начальные условия для y0 очевидны на основе закона о сохранении. Использование odeset чтобы задать опции:

  • Используйте постоянную большую матрицу, чтобы представлять левую сторону системы уравнений.

$$\left( \begin{array}{c} y'_1\\ y'_2\\ 0 \end{array} \right) = M y'
\rightarrow M = \left( \begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 &
0 \end{array} \right)$$

  • Установите допуск относительной погрешности равным 1e-4.

  • Используйте абсолютную погрешность 1e-10 для второго компонента решения, поскольку шкала существенно изменяется по сравнению с другими компонентами.

  • Выйдите из 'MassSingular' опция на своём значении по умолчанию 'maybe' для проверки автоматического обнаружения ДАУ.

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');

Входные параметры

свернуть все

Функции, которые нужно решить, заданные как указатель на функцию, который определяет функции, которые будут интегрированы.

Функция dydt = odefun(t,y), для скаляра t и вектор-столбец y, должен вернуть вектор-столбец dydt типа данных single или double который соответствует f(t,y). odefun необходимо принять оба входных параметров, t и y, даже если один из аргументов не используется в функции.

Для примера, чтобы решить y'=5y−3, используйте функцию:

function dydt = odefun(t,y)
dydt = 5*y-3;

Для системы уравнений выход odefun является вектором. Каждый элемент в векторе является решением одного уравнения. Для примера, чтобы решить

y'1=y1+2y2y'2=3y1+2y2

используйте функцию:

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);

Для получения информации о том, как предоставить дополнительные параметры функции odefun, см. Параметризация функций.

Пример: @myFcn

Типы данных: function_handle

Интервал интегрирования, заданный как вектор. Как минимум, tspan должен быть вектором с двумя элементами [t0 tf] определение начального и последнего времени. Получение решений в определенное время между t0 и tf, используйте более длинный вектор вида [t0,t1,t2,...,tf]. Элементы в tspan все должно увеличиваться или все уменьшаться.

Решатель накладывает начальные условия, заданные y0 в начальное время tspan(1), затем интегрируется из tspan(1) на tspan(end):

  • Если tspan имеет два элемента, [t0 tf]затем решатель возвращает решение, оцененное на каждом внутреннем шаге интегрирования в пределах интервала.

  • Если tspan имеет более двух элементов [t0,t1,t2,...,tf], затем решатель возвращает решение, оцененное в заданных точках. Однако решатель не шагает точно к каждой точке, заданной в tspan. Вместо этого решатель использует свои собственные внутренние шаги, чтобы вычислить решение, затем оценивает решение в запрашиваемых точках в tspan. Решения, полученные в заданных точках, имеют тот же порядок точности, что и решения, вычисленные на каждом внутреннем шаге.

    Определение нескольких промежуточных точек мало влияет на эффективность расчетов, но для больших систем это может повлиять на управление памятью.

Значения tspan используются решателем для вычисления подходящих значений для InitialStep и MaxStep:

  • Если tspan содержит несколько промежуточных точек [t0,t1,t2,...,tf], затем указанные точки дают указание на шкалу для задачи, которая может повлиять на значение InitialStep используется решателем. Поэтому решение, полученное решателем, может быть разным в зависимости от того, задаете ли вы tspan как двухэлементный вектор или как вектор с промежуточными точками.

  • Начальное и окончательное значения в tspan используются для вычисления максимального размера шага MaxStep. Поэтому изменение начальных или конечных значений в tspan может привести к тому, что решатель будет использовать другую последовательность шагов, что может изменить решение.

Пример: [1 10]

Пример: [1 3 5 7 9 10]

Типы данных: single | double

Начальные условия, заданные как вектор. y0 должна быть такой же длины, как и вектора выхода odefun, так что y0 содержит начальное условие для каждого уравнения, заданное в odefun.

Типы данных: single | double

Структура опции, заданная как массив структур. Используйте odeset функция для создания или изменения структуры опций. Список опций, совместимых с каждым решателем, см. в Сводных данных опций ОДУ.

Пример: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) задает относительную погрешность допуск 1e-5, включает отображение статистики решателя и задает выходную функцию @odeplot для построения графика решения по мере его вычисления.

Типы данных: struct

Выходные аргументы

свернуть все

Оценочные точки, возвращенные как вектор-столбец.

  • Если tspan содержит два элемента, [t0 tf], затем t содержит внутренние точки оценки, используемые для выполнения интегрирования.

  • Если tspan содержит более двух элементов, затем t то же, что и tspan.

Решения, возвращенные как массив. Каждая строка в y соответствует решению при значении, возвращенном в соответствующей строке t.

Время событий, возвращаемое как вектор-столбец. Время событий в te соответствуют решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Решение во время событий, возвращаемое как массив. Время событий в te соответствуют решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Индекс функции инициируемого события, возвращаемый как вектор-столбец. Время событий в te соответствуют решениям, возвращаемым в ye, и ie определяет, какое событие произошло.

Структура для оценки, возвращенная как массив структур. Используйте эту структуру с deval функция для вычисления решения в любой точке интервала [t0 tf]. The sol массив структур всегда включает в себя следующие поля:

Структурное полеОписание

sol.x

Вектор строка с шагом, выбранным решателем.

sol.y

Решения. Каждый столбец sol.y(:,i) содержит решение во время sol.x(i).

sol.solver

Имя решателя.

Кроме того, если вы задаете Events опция и события обнаруживаются, затем sol также включает в себя следующие поля:

Структурное полеОписание

sol.xe

Указывает, когда произошли события. sol.xe(end) содержит точную точку терминального события, если она имеется.

sol.ye

Решения, которые соответствуют событиям в sol.xe.

sol.ie

Индексы в вектор, возвращенный функцией, заданной в Events опция. Значения указывают, какое событие обнаружил решатель.

Алгоритмы

ode15s является решателем переменного шага переменного порядка (VSVO), основанным на численных формулах дифференцирования (NDF) порядков 1-5. Опционально он может использовать формулы дифференцирования назад (BDF, также известные как метод Гира), которые обычно менее эффективны. Как ode113, ode15s является многоступенчатым решателем. Использовать ode15s если ode45 отказывает или очень неэффективен, и вы подозреваете, что задача жесткая, или при решении дифференциально-алгебраического уравнения (ДАУ) [1], [2].

Ссылки

[1] шемпин, L. F. and M. W. Reichelt, «The MATLAB ODE Suite», SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1-22.

[2] шемпин, L. F., M. W. Reichelt, and J.A. Кирженка, «Solving Index-1 DAEs in MATLAB and Simulink», SIAM Review, Vol. 41, 1999, pp. 538-552.

См. также

| | | | | |

Темы

Представлено до R2006a

Документация MATLAB

Поддержка