NDSolve and hybrid dynamics (Differential Algebraic Equation DAE)
Ravi Balasubramanian
I am trying to solve a hybrid dynamics problem using NDSolve.
I am aware that I could use the event locator methods in NDSolve and
manually take care of the transitions (set initial conditions etc). But
my system is complicated, and I would have to do a lot of book-keeping.
The following represents a simplified version of the problem.
eqns = {x'[t] == y[t] + Sin[t], y[t] == If[t < 0.1, 1, 2],
x[-Pi] == 1/2};
sol = NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
It is a first order differential equation in x[t], but it includes a
y[t] variable that is discontinuous. There should be solution for this
problem, but I get the following error message from NDSolve.
LinearSolve::"nosol" : "Linear equation encountered that has no solution.
The error is probably arising from the "If" condition.
FYI, in my real system the y[t] variable is a Lagrange multiplier that
is enforcing a constraint. Under certain conditions, that constraint is
broken and the system transitions into another state. It will be great
if this problem can be solved with NDSolve keeping track of the y[t]
variable for me. Any suggestions greatly appreciated.
Thanks!
Ravi Balasubramanian
Yale University
Ravi Balasubramanian
Thanks for your suggestion.
Unfortunately, a smoothing function would not be useful, since the
system has binary states.
I figured a way to do this is to include the f[x] function into the
x'[t] equation such as
x'[t] == f[x[t]] + Sin[t]
(f[x] has only two states and no smoothing)
and use NDSolve to solve for x[t]. This works. And then I can evaluate
f[x[t]] using the NDSolve solution.
Just to summarize this discussion: NDSolve does not support binary
variables as part of the state?
Ravi
Oliver Ruebenkoenig wrote:
>
>
> Ravi,
>
> would a smoothing function be useful?
>
> f[x_] := Which[
> x <= -1, 0,
> -1 < x < 1, 1/2 + x (15/16 - x^2 (5/8 - 3/16 x^2)),
> x >= 1, 1] + 1
>
> Plot[f[(x - 0.1)/0.001], {x, -0.2, 0.2}]
>
> eqns = {x'[t] == y[t] + Sin[t], y[t] == f[t], x[-\[Pi]] == 1/2};
> sol = NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
>
> Hope this helps,
>
> Oliver
Bob Hanlon
If is a programming construct not a mathematical function. Use Piecewise
eqns = {x'[t] == y[t] + Sin[t],
y[t] == Piecewise[{{1, t < 1/10}}, 2],
x[-Pi] == 1/2};
sol = {x[t], y[t]} /. DSolve[eqns, {x[t], y[t]}, t][[1]]
{(1/2)*(-1 + 2*Pi + 2*Piecewise[{{t - Cos[t], t <= 1/10}},
-(1/10) + 2*t - Cos[t]]), Piecewise[{{1, t < 1/10}}, 2]}
Plot[sol, {t, -5, 1},
Frame -> True,
Axes -> False,
Exclusions -> 1/10,
ExclusionsStyle -> Red]
For NDSolve break the problem into two regions
Bob Hanlon
---- Ravi Balasubramanian <ravi.balas...@yale.edu> wrote:
=============
Bob Hanlon
The equations are not the same if you change y[t] to y[x[t]]. Further, if you want the break in y[t] when x[t] is 0.1 then that must be the boundary condition on x rather than x[-Pi] == 0.5
eqns = {x'[t] == y[t] + Sin[t],
y[t] == Piecewise[{{1, t < t1}}, 2],
x[t1] == 1/10};
sol = {x[t], y[t]} /.
DSolve[eqns, {x[t], y[t]}, t][[1]] //
Simplify
{1/10 - 2*t1 + Cos[t1] + Piecewise[{{t - Cos[t], t < t1}}, 2*t - Cos[t]],
Piecewise[{{1, t < t1}}, 2]}
Simplify[sol, t == t1]
{1/10, 2}
Solve[(Simplify[sol[[1]], t < t1] /. t -> t1) == 1/10, t1][[1]]
{t1 -> 0}
sol = Simplify[sol /. t1 -> 0]
{11/10 + Piecewise[{{t - Cos[t], t < 0}}, 2*t - Cos[t]],
Piecewise[{{1, t < 0}}, 2]}
Plot[sol, {t, -5, 5},
Frame -> True,
Axes -> False,
Exclusions -> 0,
ExclusionsStyle -> Red]
Bob Hanlon
---- Ravi Balasubramanian <ravi.balas...@yale.edu> wrote:
=============
Thank you, Bob. I guess I oversimplified the original problem. The
variable y is state-dependent such as
y[x] := Piecewise[{{1, x < 1/10}}, 2]
I tried the following system of equations, but it failed:
eqns = {x'[t] == y[t] + Sin[t],
y[t] == Piecewise[{{1, x[t] < 1/10}}, 2], x[-Pi] == 1/2};
sol = {x[t], y[t]} /. NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
NDSolve::ndsz: At t == -3.33251, step size is effectively zero;
singularity or stiff system suspected.
I presume this is what you meant by breaking the problem into two
regions, but that transition is dependent on x.
With your suggestion of using Piecewise, I was able to run the same
system by creating a function y
y//Clear
y[xVar_?NumericQ] := Piecewise[{{1, xVar < 1/10}}, 2]
eqns = {x'[t] == y[x[t]] + Sin[t], x[-Pi] == 1/2};
xSoln = NDSolve[eqns, {x[t]}, {t, -5, 5}][[1, 1, 2]]
Plot[xSoln, {t, -5, 5}, Frame -> True, Axes -> False,
PlotStyle -> Thick]
Plot[y[xSoln], {t, -5, 5}, Frame -> True, Axes -> False,
PlotStyle -> Thick]
Is there a way that I can get NDSolve to manage the y variable for me as
part of the state because my system will have multiple such non-smooth
variables?
Ravi
Ravi Balasubramanian
schochet123
states numerically by a smoothed function. Consider, for example the
simple continuous approximation
cpw[x_]=Piecewise[{{1, x < 1/10 - eps}, {1 + (x - (1/10 - eps))/(2
eps),
1/10 - eps <= x <= 1/10 + eps}}, 2]
Taking, for example,
eps=1/1000;
substituting this into the original DAE, i.e., solving the problem
appreqns = {x'[t] == y[t] + Sin[t], y[t] == cpw[x[t]],
x[-Pi] == 1/2};
apprsol = {x[t], y[t]} /. NDSolve[appreqns, {x[t], y[t]}, {t, -5, 5}]
yields a solution which differs from that obtained by using the true
binary state by no more than about 10^(-4). If that is not good
enough, then take a smaller eps and/or a smoother approximation.
One reason for using the smoothed function is that it is not always
easy/possible to solve the algebraic equation for the slaved variable.
Consider for example
alteqns = {x'[t] == y[t] + Sin[t],
y[t]^3 - (2 x[t] + 100) y[t] == cpw[x[t]], x[-Pi] == 1/2,
y[-Pi] == Chop[N[y /. ysubs[[1]], 50]]};
altsol = {x[t], y[t]} /. NDSolve[alteqns2, {x[t], y[t]}, {t, -5, 5}]
where
ysubs = Solve[y^3 - 101 y == 2];
In this example it is possible to omit the initial condition for y and
then NDSolve picks one of the three possible values.
Steve
On Nov 2, 12:00 pm, Ravi Balasubramanian
<ravi.balasubraman...@yale.edu> wrote:
> Hello Oliver,
>
> Thanks for your suggestion.
>
> Unfortunately, a smoothing function would not be useful, since the
> system has binary states.
>
> I figured a way to do this is to include the f[x] function into the
> x'[t] equation such as
>
> x'[t] == f[x[t]] + Sin[t]
>
> (f[x] has only two states and no smoothing)
>
> and use NDSolve to solve for x[t]. This works. And then I can evalu=
ate
> f[x[t]] using the NDSolve solution.
>
> Just to summarize this discussion: NDSolve does not support binary
> variables as part of the state?
>
> Ravi
>
> Oliver Ruebenkoenig wrote:
>
> > Ravi,
>
> > would a smoothing function be useful?
>
> > f[x_] := Which[
> > x <= -1, 0,
> > -1 < x < 1, 1/2 + x (15/16 - x^2 (5/8 - 3/16 x^2)),
> > x >= 1, 1] + 1
>
> > Plot[f[(x - 0.1)/0.001], {x, -0.2, 0.2}]
>
> > eqns = {x'[t] == y[t] + Sin[t], y[t] == f[t], x[-\[Pi]] ==
= 1/2};
> > sol = NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
>
> > Hope this helps,
>
> > Oliver
>
> > On Mon, 1 Nov 2010, Ravi Balasubramanian wrote:
>
> >> Friends,
>
> >> I am trying to solve a hybrid dynamics problem using NDSolve.
>
> >> I am aware that I could use the event locator methods in NDSolve and
> >> manually take care of the transitions (set initial conditions etc). =
But
> >> my system is complicated, and I would have to do a lot of book-keeping=
.
>
> >> The following represents a simplified version of the problem.
>
> >> eqns = {x'[t] == y[t] + Sin[t], y[t] == If[t < 0.1, 1, 2],
> >> x[-Pi] == 1/2};
> >> sol = NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
>
> >> It is a first order differential equation in x[t], but it includes a
> >> y[t] variable that is discontinuous. There should be solution for t=
his
> >> problem, but I get the following error message from NDSolve.
>
> >> LinearSolve::"nosol" : "Linear equation encountered that has no soluti=
on.
>
> >> The error is probably arising from the "If" condition.
>
> >> FYI, in my real system the y[t] variable is a Lagrange multiplier that
> >> is enforcing a constraint. Under certain conditions, that constrain=
t is
> >> broken and the system transitions into another state. It will be gr=
Albert Retey
> Thank you, Bob. I guess I oversimplified the original problem. The
> variable y is state-dependent such as
> y[x] := Piecewise[{{1, x < 1/10}}, 2]
>
> I tried the following system of equations, but it failed:
>
> eqns = {x'[t] == y[t] + Sin[t],
> y[t] == Piecewise[{{1, x[t] < 1/10}}, 2], x[-Pi] == 1/2};
> sol = {x[t], y[t]} /. NDSolve[eqns, {x[t], y[t]}, {t, -5, 5}]
>
> NDSolve::ndsz: At t == -3.33251, step size is effectively zero;
> singularity or stiff system suspected.
>
> I presume this is what you meant by breaking the problem into two
> regions, but that transition is dependent on x.
>
> With your suggestion of using Piecewise, I was able to run the same
> system by creating a function y
>
> y//Clear
> y[xVar_?NumericQ] := Piecewise[{{1, xVar < 1/10}}, 2]
> eqns = {x'[t] == y[x[t]] + Sin[t], x[-Pi] == 1/2};
>
> xSoln = NDSolve[eqns, {x[t]}, {t, -5, 5}][[1, 1, 2]]
> Plot[xSoln, {t, -5, 5}, Frame -> True, Axes -> False,
> PlotStyle -> Thick]
>
> Plot[y[xSoln], {t, -5, 5}, Frame -> True, Axes -> False,
> PlotStyle -> Thick]
>
> Is there a way that I can get NDSolve to manage the y variable for me as
> part of the state because my system will have multiple such non-smooth
> variables?
As you have seen, in principle Mathematica can handle DAEs
(Differential-Algebraic Equation), but my experience is that if you can
formulate the problem as a pure DE (or set of those), then NDSolve will
in general perform better. My guess is that it then has more Methods to
choose from and most often makes a good choice. In your case the message
about the stiff system might be a hint: while the pure DE methods will
often handle stiff systems very well (presumably by using the stiffness
switching method) the DAE methods will not be able to do that. Whether
that is a problem in principle or just caused by the way the DAE solving
methods are implemented I can't say, but for the moment the situation is
like that: NDSolve does work better for pure DEs than for DAEs in many
cases...
Fortunately, for those cases where a reformulation is possible,
Mathematica can help you do the symbolic manipulations in an automated
way by using replacement rules or whatever fits for you. Here is a very
simple example for your case:
rules = {y[x_] :> Piecewise[{{1, x < 1/10}}, 2]};
eqns = {x'[t] == y[x[t]] + Sin[t], x[-Pi] == 1/2} /. rules;
xSoln = NDSolve[eqns, {x[t]}, {t, -5, 5}][[1, 1, 2]]
hth,
albert