decic

Вычисление допустимых начальных условий для ode15i

Описание

пример

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0) использует y0 и yp0 как догадки для начальных условий полностью неявной функции odefun, содержит компоненты, заданные как fixed_y0 и fixed_yp0 как фиксированное, затем вычисляет значения для нефиксированных компонентов. Результатом является полный набор последовательных начальных условий. Новые значения yo_new и yp0_new удовлетворить odefun(t0,y0_new,yp0_new) = 0 и подходят для использования в качестве начальных условий с ode15i.

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options) также использует структуру опций options для задания значений для AbsTol и RelTol. Создайте структуру опций с помощью odeset.

[y0_new,yp0_new,resnrm] = decic(___) возвращает норму odefun(t0,y0_new,yp0_new) как resnrm. Если норма кажется чрезмерно большой, то используйте options для уменьшения относительного допуска RelTol, которое имеет значение по умолчанию 1e-3.

Примеры

свернуть все

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

0=2y1-y20=y1+y2

Эти уравнения достаточно просты, чтобы считать непротиворечивые начальные условия для переменных. Для примера, если вы исправляете y1=1, затем y2=-1 согласно второму уравнению и y1=-1/2 согласно первому уравнению. Поскольку эти значения y1, y1, и y2 удовлетворить уравнения, они непротиворечивы.

Подтвердите эти значения при помощи decic чтобы вычислить допустимые начальные условия для уравнений, фиксируя значение y1=1. Используйте догадки о y0 = [1 0] и yp0 = [0 0], которые не удовлетворяют уравнениям и, таким образом, несогласованны.

odefun = @(t,y,yp) [2*yp(1)-y(2); y(1)+y(2)];
t0 = 0;
y0 = [1 0];
yp0 = [0 0];
[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[])
y0 = 2×1

     1
    -1

yp0 = 2×1

   -0.5000
         0

Вычисление допустимых начальных условий и решение неявной ОДУ с ode15i.

Уравнение Вайссинджера

ty2(y)3-y3(y)2+t(t2+1)y-t2y=0.

Поскольку уравнение находится в общей форме f(t,y,y)=0, можно использовать ode15i функция для решения неявного дифференциального уравнения.

Кодовое уравнение

Чтобы кодировать уравнение в форме, подходящей для ode15i, вам нужно написать функцию с входами для t, y, и y который возвращает остаточное значение уравнения. Функция @weissinger кодирует это уравнение. Просмотрите файл функции.

type weissinger
function res = weissinger(t,y,yp)
%WEISSINGER  Evaluate the residual of the Weissinger implicit ODE
%
%   See also ODE15I.

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

res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

Вычисление допустимых начальных условий

The ode15i решатель требует допустимых начальных условий, то есть начальные условия, подаваемые на решатель, должны удовлетворять

f(t0,y,y)=0.

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

В этом случае исправьте начальное значение y(t0)=32 и позвольте decic вычислить допустимое начальное значение для производной y(t0), начиная с начального предположения о y(t0)=0.

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0)
y0 = 1.2247
yp0 = 0.8165

Решение уравнения

Используйте согласованные начальные условия, возвращенные decic с ode15i для решения ОДУ за временной интервал [110].

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);

Графическое изображение результатов

Точное решение этой ОДУ

y(t)=t2+12.

Постройте график численного решения y вычисляется по ode15i против аналитического решения ytrue.

ytrue = sqrt(t.^2 + 0.5);
plot(t,y,'*',t,ytrue,'-o')
legend('ode15i', 'exact')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent ode15i, exact.

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

свернуть все

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

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

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

function f = odefun(t,y,yp)
f = yp - y;

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

y'1y2=0y'2+1=0,

использовать эту функцию.

function dy = odefun(t,y,yp)
dy = zeros(2,1);
dy(1) = yp(1)-y(2);
dy(2) = yp(2)+1;

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

Пример: @myFcn

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

Начальное время, заданное как скаляр. decic использует начальное время для вычисления допустимых начальных условий, удовлетворяющих odefun(t0,y0_new,yp0_new) = 0.

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

Начальные предположения для y-компоненты, заданные как вектор. Каждый элемент в y0 задает начальное условие для одной зависимой переменной yn в системе уравнений, заданных odefun.

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

y-компоненты для удержания фиксированными, заданные как вектор на 1с и 0с, или как [].

  • Задайте fixed_y0(i) = 1 если изменение не разрешено в предположении для y0(i).

  • Задайте fixed_y0 = [] если любая запись может быть изменена.

Вы не можете исправить больше length(yp0) компоненты. В зависимости от конкретной задачи не всегда можно исправить определенные компоненты y0 или yp0. Это лучшая практика, чтобы не исправить больше компонентов, чем необходимо.

Начальные предположения для y'-компоненты, заданные как вектор. Каждый элемент в yp0 задает начальное условие для одной дифференцированной зависимой переменной y'n в системе уравнений, заданных odefun.

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

y'-компоненты для удержания фиксированными, заданные как вектор на 1с и 0с, или как [].

  • Задайте fixed_yp0(i) = 1 если изменение не разрешено в предположении для yp0(i).

  • Задайте fixed_yp0 = [] если любая запись может быть изменена.

Вы не можете исправить больше length(yp0) компоненты. В зависимости от конкретной задачи не всегда можно исправить определенные компоненты y0 или yp0. Это лучшая практика, чтобы не исправить больше компонентов, чем необходимо.

Структура опций, заданная как массив структур. Используйте odeset функция для создания или изменения структуры опций. Соответствующие опции для использования с decic функции RelTol и AbsTol, которые управляют порогами ошибок, используемыми для вычисления начальных условий.

Пример: options = odeset('RelTol',1e-5)

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

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

свернуть все

Допустимые начальные условия для y0, возвращается как вектор. Если значение resnrm является маленьким, тогда yo_new и yp0_new удовлетворить odefun(t0,y0_new,yp0_new) = 0 и подходят для использования в качестве начальных условий с ode15i.

Допустимые начальные условия для yp0, возвращается как вектор. Если значение resnrm является маленьким, тогда yo_new и yp0_new удовлетворить odefun(t0,y0_new,yp0_new) = 0 и подходят для использования в качестве начальных условий с ode15i.

Норма невязки, возвращенная как вектор. resnrm - норма odefun(t0,y0_new,yp0_new).

  • Небольшое значение resnrm указывает, что decic успешно вычисленные допустимые начальные условия, которые удовлетворяют odefun(t0,y0_new,yp0_new) = 0.

  • Если значение resnrm большой, попробуйте настроить пороги ошибок RelTol и AbsTol использование options вход.

Совет

  • The ihb1dae и iburgersode пример использования файлов decic вычисление допустимых начальных условий перед решением с ode15i. Тип edit ihb1dae или edit iburgersode чтобы просмотреть код.

  • Вы можете дополнительно использовать decic вычисление допустимых начальных условий для ДАУ, решаемых ode15s или ode23t. Для этого выполните следующие шаги.

    1. Переписать систему уравнений в полностью неявной форме f(t,y,y') = 0.

    2. Звонить decic вычислить допустимые начальные условия для уравнений.

    3. Задайте y0_new в качестве начального условия в вызове решателя и задайте yp_new как значение InitialSlope опция odeset.

См. также

| |

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