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

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

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

Найдите эти указатели на функцию при помощи 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. Поиск начальных условий

Решатели требуют начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнениям при помощи MATLAB decic функция. 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));

Найдите непротиворечивые начальные значения для переменных и их производных с помощью MATLAB decic функция. Первый аргумент 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