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

Рассмотрите 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)