Найдите последовательные начальные условия для неявной системы ОДУ первого порядка с алгебраическими ограничениями
[ находит допустимые начальные условия для системы неявных обыкновенных дифференциальных уравнений первого порядка с алгебраическими ограничениями, возвращаемыми y0,yp0]
= decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)reduceDAEToODE функция.
Вызов [eqs,constraintEqs] = reduceDAEToODE(DA_eqs,vars) уменьшает систему дифференциальных алгебраических уравнений DA_eqs в систему неявных ОДУ eqs. Это также возвращает ограничительные уравнения, встречающиеся во время сокращения системы. Для переменных этой системы ОДУ и их производных, decic находит допустимые начальные условия y0, yp0 в то время t0.
Подстановка числовых значений y0, yp0 в дифференциальные уравнения subs(eqs, [t; vars(t); diff(vars(t))], [t0; y0; yp0]) и ограничительные уравнения subs(constr, [t; vars(t); diff(vars(t))], [t0; y0; yp0]) производит нулевые векторы. Здесь, vars должен быть вектор-столбец.
y0_est определяет числовые оценки для значений переменных vars в то время t0, и fixedVars указывает значения в y0_est это не должно меняться во время численного поиска. Необязательный аргумент yp0_est позволяет задать числовые оценки для значений производных переменных vars в то время t0.
Сократите систему ДАУ до системы неявных ОДУ. Затем найдите согласованные начальные условия для переменных получившейся системы ODE и их первых производных.
Создайте следующую дифференциальную алгебраическую систему.
syms x(t) y(t)
DA_eqs = [diff(x(t),t) == cos(t) + y(t),...
x(t)^2 + y(t)^2 == 1];
vars = [x(t); y(t)];Использовать reduceDAEToODE преобразование этой системы в систему неявных ОДУ.
[eqs, constraintEqs] = reduceDAEToODE(DA_eqs, vars)
eqs =
diff(x(t), t) - y(t) - cos(t)
- 2*x(t)*diff(x(t), t) - 2*y(t)*diff(y(t), t)
constraintEqs =
1 - y(t)^2 - x(t)^2Создайте набор опций, который задает числовые допуски для численного поиска.
options = odeset('RelTol', 10.0^(-7), 'AbsTol', 10.0^(-7));Исправьте значения t0 = 0 для временных и численных оценок для допустимых значений переменных и их производных.
t0 = 0; y0_est = [0.1, 0.9]; yp0_est = [0.0, 0.0];
Можно рассматривать ограничение как алгебраическое уравнение для переменной x с фиксированным параметром y. Для этого задайте fixedVars = [0 1]. Также можно рассматривать его как алгебраическое уравнение для переменной y с фиксированным параметром x. Для этого задайте fixedVars = [1 0].
Сначала установите начальное значение x(t0) = y0_est(1) = 0.1.
fixedVars = [1 0]; [y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)
y0 =
0.1000
0.9950
yp0 =
1.9950
-0.2005Теперь измените fixedVars на [0 1]. Это исправляет y(t0) = y0_est(2) = 0.9.
fixedVars = [0 1]; [y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)
y0 =
-0.4359
0.9000
yp0 =
1.9000
0.9202Проверьте, что эти начальные значения являются допустимыми начальными значениями, удовлетворяющими уравнениям и ограничениям.
subs(eqs, [t; vars; diff(vars,t)], [t0; y0; yp0])
ans = 0 0
subs(constraintEqs, [t; vars; diff(vars,t)], [t0; y0; yp0])
ans = 0
daeFunction | findDecoupledBlocks | incidenceMatrix | isLowIndexDAE | massMatrixForm | odeFunction | reduceDAEIndex | reduceDAEToODE | reduceDifferentialOrder | reduceRedundancies