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

Шаг 1. Преобразуйте ДАУ в Указатели на функции

От вывода reduceDAEIndex у вас есть вектор уравнений DAEs и вектор переменных DAEvars. Чтобы использовать ode15s или ode23t, вам нужны два указателя на функцию: одно представление большой матрицы системы ДАУ и другого представления правых сторон уравнений большой матрицы. Эти указатели на функцию формируют эквивалентное представление большой матрицы системы ОДУ где M (t, y (t)) y’ (t) = f (t, y (t)).

Найдите, что эти указатели на функцию при помощи massMatrixForm получают большую матрицу M и правые стороны F.

[M,f] = massMatrixForm(DAEs,DAEvars)

M =
[ 0,  0, 0, 0,  0, 0, 0]
[ 0,  0, 0, 0,  0, 0, 0]
[ 0,  0, 0, 0,  0, 0, 0]
[ 0,  0, 0, 0,  0, 0, 0]
[ 0,  0, 0, 0,  0, 0, 0]
[ 0,  0, 0, 0, -1, 0, 0]
[ 0, -1, 0, 0,  0, 0, 0]

f =
                                 (T(t)*x(t) - m*r*Dxtt(t))/r
                        -(g*m*r - T(t)*y(t) + m*r*Dytt(t))/r
                                       r^2 - y(t)^2 - x(t)^2
                             - 2*Dxt(t)*x(t) - 2*Dyt(t)*y(t)
 - 2*Dxtt(t)*x(t) - 2*Dytt(t)*y(t) - 2*Dxt(t)^2 - 2*Dyt(t)^2
                                                    -Dytt(t)
                                                     -Dyt(t)

Уравнения в DAEs могут содержать символьные параметры, которые не заданы в векторе переменных DAEvars. Найдите эти параметры при помощи setdiff на выводе symvar от DAEs и DAEvars.

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars)

extraParams =
[ g, m, r]

Большая матрица M не имеет этих дополнительных параметров. Поэтому преобразуйте M непосредственно в указатель на функцию при помощи odeFunction.

M = odeFunction(M, DAEvars);

Преобразуйте f в указатель на функцию. Задайте дополнительные параметры как дополнительные входные параметры к odeFunction.

f = odeFunction(f, DAEvars, g, m, r);

Остальная часть рабочего процесса является чисто числовой. Установите значения параметров и создайте указатель на функцию.

g = 9.81;
m = 1;
r = 1;
F = @(t, Y) f(t, Y, g, m, r);

Шаг 2. Найдите начальные условия

Решатели требуют начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнения при помощи функции decic MATLAB. decic принимает предположения (который не может удовлетворить уравнения) для начальных условий, и пытается найти удовлетворительные начальные условия с помощью тех предположений. decic может перестать работать, в этом случае необходимо вручную предоставить сопоставимые начальные значения для проблемы.

Во-первых, проверяйте переменные в DAEvars.

DAEvars

DAEvars =
    x(t)
    y(t)
    T(t)
  Dxt(t)
  Dyt(t)
 Dytt(t)
 Dxtt(t)

Здесь, Dxt(t) является первой производной x(t), Dytt(t) является второй производной y(t) и так далее. Существует 7 переменных в 7-by-1 вектор. Таким образом предположениями для начальных значений переменных и их производных должен также быть 7-by-1 векторы.

Примите, что начальное угловое смещение маятника составляет 30 ° или pi/6, и источник координат в точке приостановки маятника. Учитывая, что мы использовали радиус r 1, начальное горизонтальное положение, x(t) является r*sin(pi/6). Начальным вертикальным положением y(t) является -r*cos(pi/6). Задайте эти начальные значения переменных в векторном y0est.

Произвольно установите начальные значения остающихся переменных и их производных к 0. Это не хорошие предположения. Однако они достаточны для нашей проблемы. В вашей проблеме, если ошибки decic, то обеспечивают лучшие предположения и относятся к странице decic.

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0; 0; 0];
yp0est = zeros(7,1);

Создайте набор опции, который содержит большую матрицу M и исходные предположения yp0est, и задает числовые допуски к числовому поиску.

opt = odeset('Mass', M, 'InitialSlope', yp0est,...
                'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));

Найдите сопоставимые начальные значения для переменных и их производных при помощи функции decic MATLAB. Первый аргумент decic должен быть указателем на функцию, описывающим ДАУ как f(t,y,yp) = f(t,y,y') = 0. С точки зрения M и F, это означает f(t,y,yp) = M(t,y)*yp - F(t,y).

implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
[y0, yp0] = decic(implicitDAE, 0, y0est, [], yp0est, [], opt)

y0 =
    0.4771
   -0.8788
   -8.6214
         0
    0.0000
   -2.2333
   -4.1135

yp0 =
         0
    0.0000
         0
         0
   -2.2333
         0
         0

Теперь создайте набор опции, который содержит большую матрицу M системы и векторный yp0 сопоставимых начальных значений для производных. Вы будете использовать этот набор опции при решении системы.

opt = odeset(opt, 'InitialSlope', yp0);

Шаг 3. Решите систему ДАУ

Решите систему, объединяющуюся по отрезку времени 0 ≤ t ≤ 0.5. Добавьте линии сетки и легенду к графику. Код использует ode15s. Вместо этого можно использовать ode23s, заменяя ode15s на ode23s.

[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt);
plot(tSol,ySol(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(DAEvars(k));
end

legend(S, 'Location', 'Best')
grid on

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

Установите r на 2 и регенерируйте указатель на функцию и начальные условия.

r = 2;

F = @(t, Y) f(t, Y, g, m, r);
y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0; 0; 0];
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
[y0, yp0] = decic(implicitDAE, 0, y0est, [], yp0est, [], opt);

opt = odeset(opt, 'InitialSlope', yp0);

Решите систему для нового значения параметров.

[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt);
plot(tSol,ySol(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(DAEvars(k));
end

legend(S, 'Location', 'Best')
grid on