Решите ДАУ Используя решатели большой матрицы

Решите дифференциальные алгебраические уравнения при помощи одного из решателей большой матрицы, доступных в MATLAB®. Чтобы использовать этот рабочий процесс, сначала полные шаги 1, 2, и 3 от Решают Дифференциальные Алгебраические уравнения (ДАУ). Затем используйте решатель большой матрицы вместо ode15i.

Этот пример демонстрирует использование ode15s или ode23t. Для получения дополнительной информации на других решателях, смотрите, Выбирают ODE Solver (MATLAB) и адаптируют рабочий процесс на этой странице.

Шаг 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 переменных- 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));

Найдите сопоставимые начальные значения для переменных и их производных при помощи 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. Решите систему ДАУ

Решите систему, объединяющуюся по отрезку времени 0t≤ 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

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте