Find consistent initial conditions for first-order implicit ODE system with algebraic constraints
[
finds
consistent initial conditions for the system of first-order implicit
ordinary differential equations with algebraic constraints returned
by the y0
,yp0
]
= decic(eqs
,vars
,constraintEqs
,t0
,y0_est
,fixedVars
,yp0_est
,options
)reduceDAEToODE
function.
The call [eqs,constraintEqs] = reduceDAEToODE(DA_eqs,vars)
reduces
the system of differential algebraic equations DA_eqs
to
the system of implicit ODEs eqs
. It also returns
constraint equations encountered during system reduction. For the
variables of this ODE system and their derivatives, decic
finds
consistent initial conditions y0
, yp0
at
the time t0
.
Substituting the numerical values y0
, yp0
into
the differential equations subs(eqs, [t; vars(t); diff(vars(t))],
[t0; y0; yp0])
and the constraint equations subs(constr,
[t; vars(t); diff(vars(t))], [t0; y0; yp0])
produces zero
vectors. Here, vars
must be a column vector.
y0_est
specifies numerical estimates for
the values of the variables vars
at the time t0
,
and fixedVars
indicates the values in y0_est
that
must not change during the numerical search. The optional argument yp0_est
lets
you specify numerical estimates for the values of the derivatives
of the variables vars
at the time t0
.
Reduce the DAE system to a system of implicit ODEs. Then, find consistent initial conditions for the variables of the resulting ODE system and their first derivatives.
Create the following differential algebraic system.
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)];
Use reduceDAEToODE
to convert this system
to a system of implicit ODEs.
[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
Create an option set that specifies numerical tolerances for the numerical search.
options = odeset('RelTol', 10.0^(-7), 'AbsTol', 10.0^(-7));
Fix values t0 = 0
for the time and numerical
estimates for consistent values of the variables and their derivatives.
t0 = 0; y0_est = [0.1, 0.9]; yp0_est = [0.0, 0.0];
You can treat the constraint as an algebraic equation for the
variable x
with the fixed parameter y
.
For this, set fixedVars = [0 1]
. Alternatively,
you can treat it as an algebraic equation for the variable y
with
the fixed parameter x
. For this, set fixedVars
= [1 0]
.
First, set the initial value 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
Now, change fixedVars
to [0 1]
.
This fixes 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
Verify that these initial values are consistent initial values satisfying the equations and the constraints.
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