Решите дифференциальные алгебраические уравнения при помощи одного из решателей большой матрицы, доступных в MATLAB®. Чтобы использовать этот рабочий процесс, сначала полные шаги 1, 2, и 3 от Решают Дифференциальные Алгебраические уравнения (ДАУ). Затем используйте решатель большой матрицы вместо ode15i
.
Этот пример демонстрирует использование ode15s
или ode23t
. Для получения дополнительной информации на других решателях, смотрите, Выбирают ODE Solver (MATLAB) и адаптируют рабочий процесс на этой странице.
От выхода 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);
Решатели требуют начальных значений для всех переменных в указателе на функцию. Найдите начальные значения, которые удовлетворяют уравнениям при помощи 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);
Решите систему, объединяющуюся по отрезку времени 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