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

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

Рассмотрите 2D физический маятник, состоя из массового 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 = 
$$\left(\begin{array}{c}m\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }x\left(t\right)=\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ m\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }y\left(t\right)=\frac{F\left(t\right)\hspace{0.17em}y\left(t\right)}{r}-g\hspace{0.17em}m\\ {x\left(t\right)}^2+{y\left(t\right)}^2=r^2\end{array}\right)$$

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

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

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

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

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

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

newVars = 
$$\left(\begin{array}{cc}\mathrm{Dxt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ \mathrm{Dyt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

Попытайтесь решить систему ДАУ Высокого Индекса

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

Чтобы преобразовать систему ДАУ в указатель функции 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, чтобы установить числовые допуски. Затем используйте функцию decic MATLAB, чтобы вычислить сопоставимые начальные условия 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);

warning(s)

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

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

isLowIndexDAE(eqs, vars)

ans = logical
   0

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

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

eqs = 
$$\left(\begin{array}{c}m\hspace{0.17em}\mathrm{Dxtt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ g\hspace{0.17em}m+m\hspace{0.17em}\mathrm{Dytt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}y\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ \mathrm{Dxt}\left(t\right)-{\mathrm{Dxt}}_1\left(t\right)\\ \mathrm{Dyt}\left(t\right)-{\mathrm{Dyt}}_1\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}}_1\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}{\mathrm{Dyt}}_1\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}y\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)+2\hspace{0.17em}{{\mathrm{Dxt}}_1\left(t\right)}^2+2\hspace{0.17em}{{\mathrm{Dyt}}_1\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxt1t}\left(t\right)\hspace{0.17em}x\left(t\right)\\ \mathrm{Dxtt}\left(t\right)-\mathrm{Dxt1t}\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

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

newVars = 
$$\left(\begin{array}{cc}\mathrm{Dytt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }\mathrm{Dxt}\left(t\right)\\ {\mathrm{Dxt}}_1\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\\ \mathrm{Dxt1t}\left(t\right)& \frac{{\partial }^2}{\partial t^2}\mathrm{ }x\left(t\right)\end{array}\right)$$

index = 3

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

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

eqs = 
$$\left(\begin{array}{c}-\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)-m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxtt}\left(t\right)}{r}\\ \frac{g\hspace{0.17em}m\hspace{0.17em}r-F\left(t\right)\hspace{0.17em}y\left(t\right)+m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dytt}\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ 2\hspace{0.17em}\mathrm{Dxt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dyt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2+2\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxtt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dytt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

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

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

Вычислите сопоставимые начальные условия для индекса, уменьшаемого функцией decic MATLAB. Здесь, 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)