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

Шаг 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. Найти начальные условия

Для решателей требуются начальные значения для всех переменных в дескрипторе функции. Поиск начальных значений, удовлетворяющих уравнениям, с помощью 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около-1 вектор. Таким образом, предположения для начальных значений переменных и их производных также должны быть 7около-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