Анализируйте и манипулируйте дифференциальными алгебраическими уравнениями

Открыть Live Script

Этот пример показывает, как решить дифференциальные алгебраические уравнения (ДАУ) высокого дифференциального индекса с помощью Symbolic Math Toolbox™.

Инженеры часто определяют поведение своих физических объектов (механических систем, электрических устройств и так далее) смесью дифференциальных уравнений и алгебраических уравнений. MATLAB ® предоставляет специальные числовые решатели, такие как ode15i и ode15s, способные интегрировать такие ДАУ - при условии, что их «дифференциальный индекс» не превышает 1.

Этот пример показывает рабочий процесс от настройки модели как системы дифференциальных уравнений с алгебраическими ограничениями до численной симуляции. Используются следующие функции Symbolic Math Toolbox.

  • daeFunction

  • findDecoupledBlocks

  • incidenceMatrix

  • isOfLowDAEIndex

  • reduceDifferentialOrder

  • massMatrixForm

  • reduceDAEIndex

  • reduceDAEToODE

  • reduceRedundancies

  • sym/decic

Задайте параметры модели

Рассмотрим 2-D физический маятник, состоящий из массы m присоединен к источнику строкой постоянной длины r. Только ускорение свободного падения g = 9.81 m/s^2 воздействует на массу. Модель состоит из дифференциального уравнения второго порядка для положения (x(t), y(t)) массы с неизвестной силой F(t) внутри строки, которая служит для удержания массы на круге. Сила направлена вдоль строки.

syms x(t) y(t) F(t) m g r
eqs = [m*diff(x(t), t, t) == F(t)/r*x(t);
       m*diff(y(t), t, t) == F(t)/r*y(t) - m*g;
       x(t)^2 + y(t)^2 == r^2]
eqs = 

(m ∂2∂t2 x(t)=F(t) x(t)rm ∂2∂t2 y(t)=F(t) y(t)r-g mx(t)2+y(t)2=r2)[m * diff (x (t), t, 2) = (F (t) * x (t) )/r; m * diff (y (t), t, 2) = (F (t) * y (t) )/r - g * m; x (t) ^ 2 + y (t) ^ 2 = = r ^ 2]
vars = [x(t), y(t), F(t)]
vars = (x(t)y(t)F(t))[x (t), y (t), F (t)]

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

[eqs, vars, newVars] = reduceDifferentialOrder(eqs, vars)
eqs = 

(m ∂∂t Dxt(t)-F(t) x(t)rg m+m ∂∂t Dyt(t)-F(t) y(t)r-r2+x(t)2+y(t)2Dxt(t)-∂∂t x(t)Dyt(t)-∂∂t y(t))[m * diff (Dxt (t), t) - (F (t) * x (t) )/r; g * m + m * diff (Dyt (t), t) - (F (t) * y (t) )/r; - r ^ 2 + x (t) ^ 2 + y (t) ^ 2; Dxt (t) - diff (x (t), t); Dyt (t) - diff (y (t), t)]
vars = 

(x(t)y(t)F(t)Dxt(t)Dyt(t))[x (t); y (t); F (t); Dxt (t); Краситель (t)]
newVars = 

(Dxt(t)∂∂t x(t)Dyt(t)∂∂t y(t))[Dxt (t), diff (x (t), t); Dyt (t), diff (y (t), t)]

Попробуйте решить систему ДАУ с высоким индексом

Прежде чем вы сможете использовать численный решатель MATLAB, такой как ode15i, вы должны следовать этим шагам.

  1. Преобразуйте систему ДАУ в указатель на функцию MATLAB.

  2. Выберите числовые значения для символьных параметров системы.

  3. Установите непротиворечивые начальные условия.

Чтобы преобразовать систему ДАУ в указатель на функцию MATLAB, используйте daeFunction.

F = daeFunction(eqs, vars, [m, g, r])
F = function_handle with value:
    @(t,in2,in3,in4)[in3(4,:).*in4(:,1)-(in2(3,:).*in2(1,:))./in4(:,3);in3(5,:).*in4(:,1)+in4(:,1).*in4(:,2)-(in2(3,:).*in2(2,:))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:)-in3(1,:);in2(5,:)-in3(2,:)]

Присвойте числовые значения символьным параметрам системы: m = 1kg, g = 9.18m/s^2, и r = 1m.

f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])
f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

Указатель на функцию f является подходящим входом для численного решателя ode15i. Следующим шагом является вычисление допустимых начальных условий. Использование odeset для установки числовых допусков. Затем используйте MATLAB decic функция для вычисления допустимых начальных условий y0, yp0 для позиций и производных в момент времени t0 = 0.

opt = odeset('RelTol', 10.0^(-4), 'AbsTol' , 10.0^(-4));
t0 = 0;
[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(3,1)], [], zeros(5,1), [], opt)
y0 = 5×1

    0.9777
   -0.2100
         0
         0
         0


yp0 = 5×1

         0
         0
         0
         0
   -9.8100

Протестируйте начальные условия:

f(t0, y0, yp0)
ans = 5×1
10-16 ×

         0
         0
   -0.3469
         0
         0

Теперь можно использовать ode15i чтобы попытаться решить систему. Когда вы звоните ode15i, интегрирование немедленно останавливается и выдает следующие предупреждения.

Warning: Matrix is singular, close to singular or badly scaled.

Results may be inaccurate. RCOND = NaN.

Warning: Failure at t=0.000000e+00.

Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (0.000000e+00) at time t.

В данном примере ode15i выдает эти предупреждения несколько раз. Для удобства чтения отключите предупреждения при помощи warning('off','all') перед вызовом ode15i и включите их снова.

tfinal = 0.5;
s = warning('off','all');
ode15i(f, [t0, tfinal], y0, yp0, opt);

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

warning(s)

Анализируйте и корректируйте систему ДАУ

Проверьте дифференциальный индекс системы ДАУ.

isLowIndexDAE(eqs, vars)
ans = logical
   0

Этот результат объясняет, почему ode15i невозможно решить эту систему. Эта функция требует, чтобы входная система ДАУ имела дифференциальный индекс 0 или 1. Уменьшите дифференциальный индекс путем расширения модели до эквивалентной большей системы ДАУ, которая включает некоторые скрытые алгебраические ограничения.

[eqs, vars, newVars, index] = reduceDAEIndex(eqs, vars)
eqs = 

(m Dxtt(t)-F(t) x(t)rg m+m Dytt(t)-F(t) y(t)r-r2+x(t)2+y(t)2Dxt(t)-Dxt1(t)Dyt(t)-Dyt1(t)2 Dxt1(t) x(t)+2 Dyt1(t) y(t)2 y(t) ∂∂t Dyt1(t)+2 Dxt1(t)2+2 Dyt1(t)2+2 Dxt1t(t) x(t)Dxtt(t)-Dxt1t(t)Dytt(t)-∂∂t Dyt1(t)Dyt1(t)-∂∂t y(t))[m * Dxtt (t) - (F (t) * x (t) )/r; g * m + m * Dytt (t) - (F (t) * y (t) )/r; - r ^ 2 + x (t) ^ 2 + y (t) ^ 2; Dxt (t) - Dxt1 (t); Dyt (t) - Dyt1 (t); 2 * Dxt1 (t) * x (t) + 2 * Dyt1 (t) * y (t); 2 * y (t) * diff (Dyt1 (t), t) + 2 * Dxt1 (t) ^ 2 + 2 * Dyt1 (t) ^ 2 + 2 * Dxt1t (t) * x (t); Dxtt (t) - Dxt1t (t); Dytt (t) - diff (Dyt1 (t), t); Dyt1 (t) - diff (y (t), t)]
vars = 

(x(t)y(t)F(t)Dxt(t)Dyt(t)Dytt(t)Dxtt(t)Dxt1(t)Dyt1(t)Dxt1t(t))[x (t); y (t); F (t); Dxt (t); Dyt (t); Dytt (t); Dxtt (t); Dxt1 (t); Dyt1 (t); Dxt1t (t)]
newVars = 

(Dytt(t)∂∂t Dyt(t)Dxtt(t)∂∂t Dxt(t)Dxt1(t)∂∂t x(t)Dyt1(t)∂∂t y(t)Dxt1t(t)∂2∂t2 x(t))[Dytt (t), diff (Dyt (t), t); Dxtt (t), diff (Dxt (t), t); Dxt1 (t), diff (x (t), t); Dyt1 (t), diff (y (t), t); Dxt1t (t), diff (x (t), t, 2)]
index = 3

Четвертый выход показывает, что дифференциальный индекс исходной модели равен трем. Упростите новую систему.

[eqs, vars, S] = reduceRedundancies(eqs, vars)
eqs = 

(-F(t) x(t)-m r Dxtt(t)rg m r-F(t) y(t)+m r Dytt(t)r-r2+x(t)2+y(t)22 Dxt(t) x(t)+2 Dyt(t) y(t)2 Dxt(t)2+2 Dyt(t)2+2 Dxtt(t) x(t)+2 Dytt(t) y(t)Dytt(t)-∂∂t Dyt(t)Dyt(t)-∂∂t y(t))[- (F (t) * x (t) - m * r * Dxtt (t) )/r; (g * m * r - F (t) * y (t) + m * r * Dytt (t) )/r; - r ^ 2 + x (t) ^ 2 + y (t) ^ 2; 2 * Dxt (t) * x (t) + 2 * Dyt (t) * y (t); 2 * Dxt (t) ^ 2 + 2 * Dyt (t) ^ 2 + 2 * Dxtt (t) * x (t) + 2 * Dytt (t) * y (t); Dytt (t) - diff (Dyt (t), t); Dyt (t) - diff (y (t), t)]
vars = 

(x(t)y(t)F(t)Dxt(t)Dyt(t)Dytt(t)Dxtt(t))[x (t); y (t); F (t); Dxt (t); Dyt (t); Dytt (t); Dxtt (t)]
S = struct with fields:
      solvedEquations: [3x1 sym]
    constantVariables: [0x2 sym]
    replacedVariables: [3x2 sym]
       otherEquations: [0x1 sym]

Проверьте, имеет ли новая система низкий дифференциальный индекс (0 или 1).

isLowIndexDAE(eqs, vars)
ans = logical
   1

Решение системы ДАУ с низким индексом

Сгенерируйте указатель на функцию MATLAB, который заменяет символьные параметры числовыми значениями.

F = daeFunction(eqs, vars, [m, g, r]) 
F = function_handle with value:
    @(t,in2,in3,in4)[-(in2(3,:).*in2(1,:)-in2(7,:).*in4(:,1).*in4(:,3))./in4(:,3);(-in2(3,:).*in2(2,:)+in2(6,:).*in4(:,1).*in4(:,3)+in4(:,1).*in4(:,2).*in4(:,3))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:).*in2(1,:).*2.0+in2(5,:).*in2(2,:).*2.0;in2(7,:).*in2(1,:).*2.0+in2(6,:).*in2(2,:).*2.0+in2(4,:).^2.*2.0+in2(5,:).^2.*2.0;in2(6,:)-in3(5,:);in2(5,:)-in3(2,:)]


f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])
f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

Вычисление допустимых начальных условий для индекса, уменьшенного MATLAB decic функция. Здесь, opt - структура опций, устанавливающая числовые допуски. Вы уже вычисляли его, используя odeset.

[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(5,1)], [], zeros(7,1), [], opt)
y0 = 7×1

    0.9779
   -0.2093
   -2.0528
   -0.0000
         0
   -9.3804
   -2.0074


yp0 = 7×1

         0
         0
         0
         0
   -9.3804
         0
         0

Решить систему и построить график решения.

ode15i(f, [t0, tfinal], y0, yp0, opt)

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

Документация Symbolic Math Toolbox

Поддержка

Для просмотра документации необходимо авторизоваться на сайте