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

Открыть сценарий в реальном времени

В этом примере показано, как решить дифференциальные алгебраические уравнения (DAE) высокого дифференциального индекса с помощью символьных математических Toolbox™.

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

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

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 = 
 $(\begin{array}{c} m \frac{\partial^{2}}{\partial t^{2}} x (t) = \frac{F (t) x (t)}{r} \\ m \frac{\partial^{2}}{\partial t^{2}} y (t) = \frac{F (t) y (t)}{r} - g m \\ {x (t)}^{2} + {y (t)}^{2} = r^{2} \end{array})$

vars = [x(t), y(t), F(t)]

vars =  $(\begin{array}{c} x (t) & y (t) & F (t) \end{array})$

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

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

eqs = 
 $(\begin{array}{c} m \frac{\partial}{\partial t} Dxt (t) - \frac{F (t) x (t)}{r} \\ g m + m \frac{\partial}{\partial t} Dyt (t) - \frac{F (t) y (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ Dxt (t) - \frac{\partial}{\partial t} x (t) \\ Dyt (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \end{array})$

newVars = 
 $(\begin{array}{c} Dxt (t) & \frac{\partial}{\partial t} x (t) \\ Dyt (t) & \frac{\partial}{\partial t} y (t) \end{array})$

Попробуйте решить проблему с системой High-Index DAE

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

Преобразуйте систему дисковых полок в дескриптор функции MATLAB.
Выберите числовые значения для символьных параметров системы.
Установите непротиворечивые начальные условия.

Чтобы преобразовать систему дисковых полок в дескриптор функции 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)

Проверьте исходные условия:

f(t0, y0, yp0)

Теперь вы можете использовать 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 = 
 $(\begin{array}{c} m Dxtt (t) - \frac{F (t) x (t)}{r} \\ g m + m Dytt (t) - \frac{F (t) y (t)}{r} \\ - r^{2} + {x (t)}^{2} + {y (t)}^{2} \\ Dxt (t) - {Dxt}_{1} (t) \\ Dyt (t) - {Dyt}_{1} (t) \\ 2 {Dxt}_{1} (t) x (t) + 2 {Dyt}_{1} (t) y (t) \\ 2 y (t) \frac{\partial}{\partial t} {Dyt}_{1} (t) + 2 {{Dxt}_{1} (t)}^{2} + 2 {{Dyt}_{1} (t)}^{2} + 2 Dxt1t (t) x (t) \\ Dxtt (t) - Dxt1t (t) \\ Dytt (t) - \frac{\partial}{\partial t} {Dyt}_{1} (t) \\ {Dyt}_{1} (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \\ {Dxt}_{1} (t) \\ {Dyt}_{1} (t) \\ Dxt1t (t) \end{array})$

newVars = 
 $(\begin{array}{c} Dytt (t) & \frac{\partial}{\partial t} Dyt (t) \\ Dxtt (t) & \frac{\partial}{\partial t} Dxt (t) \\ {Dxt}_{1} (t) & \frac{\partial}{\partial t} x (t) \\ {Dyt}_{1} (t) & \frac{\partial}{\partial t} y (t) \\ Dxt1t (t) & \frac{\partial^{2}}{\partial t^{2}} x (t) \end{array})$

index = 3

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

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

eqs = 
 $(\begin{array}{c} - \frac{F (t) x (t) - m r Dxtt (t)}{r} \\ \frac{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) - \frac{\partial}{\partial t} Dyt (t) \\ Dyt (t) - \frac{\partial}{\partial t} y (t) \end{array})$

vars = 
 $(\begin{array}{c} x (t) \\ y (t) \\ F (t) \\ Dxt (t) \\ Dyt (t) \\ Dytt (t) \\ Dxtt (t) \end{array})$

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)

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

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

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

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

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

Определение параметров модели

Попробуйте решить проблему с системой High-Index DAE

Анализ и настройка системы дисковых полок

Решение системы дисковых полок с низким индексом

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

Поддержка