Are there any algorithms out there for solving differential equations where
the derivatives are discontinuous? I am attempting to simulate systems that

include friction and backlash. The problem is that the models for both
friction and backlash are, for practical purposes, both infinitely stiff and
with variable structure (in both cases there is at least one state that is
sometimes algebraically dependant on another and sometimes not). While I am
finding some great ways to crash MatLab and Octave, I'm not managing to
simulate my systems in a really satisfactory manner.
I'm thinking that their ought to be a solver out there that understands such
infinitely stiff systems, and which will solve right up to the
discontinuity, then solve on the "other side" of the discontinuity, perhaps
with another system discription (to account for the variable structure),
without ever trying to resolve the behavior of the system _at_ the
discontinuity.
I know from experience that if I'm thinking it's a clever idea _now_ that
it's already been done, probably before 1950. Does this ring a bell with
anyone? Do you have some names of some methods that I might use, or books
or articles that I might read?
Thanks in advance.
--------------------------------------
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
begin 666 Tim Wescott.vcf
M0D5'24XZ5D-!4D0-"E9%4E-)3TXZ,BXQ#0I..E=E<V-O='0[5&EM#0I&3CI4
M:6T@5V5S8V]T= T*14U!24P[4%)%1CM)3E1%4DY%5#IT:6U =V5S8V]T=&1E
K<VEG;BYC;VT-"E)%5CHR,# T,#,P-%0Q-S4T,#!:#0I%3D0Z5D-!4D0-"@``
`
end

The makers of gPROMS (http://www.psenterprise.com /) claim that their
software handles discontinuities. In my experience they can slow
things down, but I guess that's inevitable.
So there's at least one algorithm :)
Also Simulink should handle discontinuities. See the "Using Simulink"
document at http://www.mathworks.com /.
L.

where
Hmm. What they say contradicts my experience, and I'm not a current owner
of MatLab (I'm waiting for a customer to come along with a project to
justify it). I know that when I've tried this before, with blocks that
should handle it, Simulink just stops at the discontinuity and that's that.
I'll have to try it again when I have Matlab.
Thanks.

Tim,
That is my experience with Simulink, too. It just stops dead.
You are right on the money on the way to handle the discontinuities. There
are two commercial codes that implement these methods that I use: AMESim
("discontinuity handling") and EASY5 ("switch states"). They interpolate to
solve for the intersection and then do the bounce, so to speak, switching
the velocity direction rather than shrinking the step size down to nothing
and turning the trajectory around in a real tight continuous curve. I like
the switch states terminology personally. The infinitely stiff concept is
discussed in the classic text (it says so on the back!) on DAE -
differential algebraic equations. It is the opposite limiting behavior from
what you are asking about -- discontinuities. You need DASSL or DASPK or
RADAU and the switch states to go fast through the bumps. Simulink doesn't
have the switch states, so it goes down to machine precision time steps and
unless you're God ("The mountains flowed before the Lord" - the Prophetess
Deborah) it seems like a long time to wait.
William Gear and Linda Petzold are the names I remember off hand for
searching. Search on all the all caps words, for that matter. I think you
can download white papers on the switch states or discontinuity handling
from the respective websites for more details (and corrections, I am sure).
Todd
writes:

where
that
and
am
such
perhaps
I don't have Matlab and Mathcad does not handle limits so I roll my own. I
use a system of differential equations with rules to check to make sure the
changes do not go out of bounds. I also have a lot of 'if' statements to
handle things like static friction and dynamic friction. Also pressures
can't go below 0 and positions can't exceed the length of the cylinder..
I use runge kutta to do the actual iterations. It takes only a few pages
to do a valve and a cylinder class for hydraulic simulations.
Peter Nachtwey

>Are there any algorithms out there for solving differential equations where
>the derivatives are discontinuous? I am attempting to simulate systems that
>include friction and backlash. The problem is that the models for both
>friction and backlash are, for practical purposes, both infinitely stiff and
>with variable structure (in both cases there is at least one state that is
>sometimes algebraically dependant on another and sometimes not). While I am
>finding some great ways to crash MatLab and Octave, I'm not managing to
>simulate my systems in a really satisfactory manner.
>
>I'm thinking that their ought to be a solver out there that understands such
>infinitely stiff systems, and which will solve right up to the
>discontinuity, then solve on the "other side" of the discontinuity, perhaps
>with another system discription (to account for the variable structure),
>without ever trying to resolve the behavior of the system _at_ the
>discontinuity.
>
>I know from experience that if I'm thinking it's a clever idea _now_ that
>it's already been done, probably before 1950. Does this ring a bell with
>anyone? Do you have some names of some methods that I might use, or books
>or articles that I might read?
>
>Thanks in advance.
>
the usual way to deal with discontinuities is to detect them using some
indicator functions (e.g. sign changes of some functions), try to locate
the time point of this precisely, stop integration there and immediately
resume taking the current solution as the new initial value.
the NAG library has such a code.
MATLAB's odesuite has the "event option" to do that.
otherwise.. design one yourself. this is not too hard . if you can characterize
these points by functions changing their sign there, then compute these functions
along with the solution and try "to look ahead' for the zero.
inverse interpolation using the stored function values is a good means
to detect those coming zeroes ahaed (e.g. with inverse quadratic and cubic
interpolation I once had much success ).
you must of course use an integrator with stepsize control.
a strong discontinuity will reveal itself by a very strong reduction of the
stepsize.
a very crude alternative solution to deal with the situation is to provide
a minimal stepsize and step over the time axis at least with this one.
this will jump over the point of discontinuity, producing some error in the order
of this step and lateron the stepsize will grow up automatically.
hth
peter

where
that
and
am
such
perhaps
Matlab can definitely do this with the ode suite (at least using ode45). In
your system equations you can have conditionals (if statements) which allow
for discrete system input. You can reference time or any system variables
in those conditionals, and plug the output directly into a differential.
tim

Don't try to automate this with a general solver. You'll just get
frustrated.
Lots of physics books solve ode's with discontinuous coefficients or
even delta-functions in the coefficients. Basically what you do is
integrate by hand in a small region around the bad point, then use
boundary conditions to join the solutions. For example,
y'' + ( k*2 - A*\delta (x) ) * y = 0
(here the bad point is x=0). Integrate from -a to a (a is small) and get
y'(a) - y'(-a) + 2*a*k^2 *y(0) = A * y(0)
so the first derivative is discontinuous, and the discontinuity is known,
in the limit a -> 0. Since the discontinuity is a finite jump, y(x) is
continuous at x=0.
In this simple case we can solve by hand: let
y(x<0) = r*cos(kx) + s*sin(kx)
y(x>0) = t*cos(kx) + u*sin(kx) ;
these clearly solve the equation on either side of x=0. Now we have
r = t ( y(0-) = y(0+) )
k*(u-s) = A*r
Hence the solution is
r*cos(kx) + s * sin(kx) , x < 0
r*cos(kx) + (s+A*r/k) * sin(kx) , x > 0 .
Clearly r and s are arbitrary and must be obtained from initial conditions
or boundary conditions, etc. (That is, a 2nd order equation has 2 arbitrary
constants.)
For more general cases you will have to actually integrate up to the bad
spots to get the jump conditions, then continue integrating on the
other side. You can't really expect a canned algorithm to know where
to join the solutions, so you will have to cook it up yourself. Let
me somewhat vaingloriously recommend my own lecture notes on ode's
as a starting point: look on the page
http://www.phys.virginia.edu/classes/551.jvn.fall01 /
look at the lecture notes page, and download the ode notes.

I was hoping that there was some general-purpose way of doing this,
basically an algorithm that lets you tell it that it's crossed a
discontinuity, or approximately where the discontinuity is in state space so
it can integrate up to it then tunnel over to the other side. Since I'm
expecting to be doing this sort of thing a _lot_ I can't use hand methods.
It sounds from other posts like I may have to spring for Simulink (sigh).

Look here:
http://portal.acm.org/citation.cfm?idR64&dl ¬M&coll=GUIDE
A robust procedure for discontinuity handling in continuous system
simulation, Transactions of the Society for Computer Simulation
International, Volume 2 , Issue 3 (September 1985)
L. G. Birta, T. I. Oren, D. L. Kettenis
There are more papers published after 1985, published (as I recall)
in\ ACM Transactions on Mathematical Software, ACM Transactions on
Modeling and Computer Simulation, Artificial Intelligence, and
Simulation. Unfortunately, I cannot find references at this moment.
The issue is briefly discussed in the book "Continuous System
Modeling" By Cellier, section 2.8, "Discontinuity handling".
The most up-to-date presentation of continuous-discrete system
simulation is in the book "Theory of Modeling and Simulation:
Integrating Discrete Event and Continuous Complex Dynamic Systems"
by Zeigler, Praehofer and Kim. I recommend also Praehofer's Ph.D.
thesis. It augments the book. However, the dissertation ("System
theoretic foundations for combined discrete-continuous system
simulation") is not available on line, you have to ask Dr. Praehofer
http://www.cast.uni-linz.ac.at/People/AcademicStaff/hp/Praehofer.htm#Aktuelles
On historical note:
Complete discrete-continuous simulation system is described in the
book "The GASP IV Simulation Language" by Alan Pritsker.
Discontinuities are handled by using the "event function" and
locating discontinuities by solving algebraic equation.
Unfortunately, all is in Fortran IV, and publication date is 1974.
The book is old, but the methodology described there is still in
use.
Detailed algorithm for modeling dynamical systems with
discontinuities is presented in even older book "Computer Simulation
of Dynamical Systems" by Kochenburger, chapter 6, Simulation of
Discontinuous Relations. Publication date is 1972. Half of the book
is about analog computers, but detailed algorithms for computer
simulation are also there.
A.L.

Why not use the Laplace transform instead? It is, of course, equivalent
but a more straightforward way to obtain the solution.

That should be k^2.
Using Mathematica, we take the Laplace transform of the differential
equation (generalized to include a shifted delta function)
LaplaceTransform[y''[x]+(k^2 - A DiracDelta[x-a]) y[x] == 0, x, s] //
Simplify
to get
(k^2+s^2) LaplaceTransform[y[x], x, s] == s y[0]+A y[a] E^(-a s)+y'[0]
Then algebraically solve (which is trivial),
Solve[%, LaplaceTransform[y[x], x, s]];
and take the inverse transform (this is the difficult step for more
complicated differential equations -- there are, however, arbitrary
precision methods useful for computing numerical inverses),
InverseLaplaceTransform[LaplaceTransform[y[x],x,s] /. First[%], s, x];
Putting a->0 and collecting terms
Collect[% /. a -> 0, {Cos[_], Sin[_]}, Simplify]
we get
Cos[k x] y[0] + Sin[k x] (A UnitStep[x] y[0] + y'[0])/k
which is equivalent to the solution you give.
Cheers,
Paul

--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014

Greetings All,
At this point, it may be good to discuss discontinuities in continuous
behavior in a little more detail as there are different aspects addressed
throughout this thread.
Typically the occurrences of discontinuous changes are modeled by
'zero-crossing functions'. When such a function changes its sign, some
discontinuity may be effected. For numerical solvers, this is important
information, as they tend to rely on continuity assumptions, and so any
discontinuous behavior (even if the discontinuity occurs in a higher order
derivative with respect to time) may degrade their performance.
The following needs to be considered (I will refrain from going into the
details, they can be found along with references in
http://www.robotic.dlr.de/control/staff/pjm/papers/hs99a/p.html ):
1. Event detection: The sign change of the zero-crossing function is
evaluated at integration points. If the function has an even number of zeros
in the interval between two integration points, you would miss it.
2. Event location: Once you know a sign change occurs between two points in
time, you may want to locate exactly where the zero crossing occurs (note
that there are applications where you rather not do this).
After these two activities, you can safely make the change to your model and
restart your numerical integration. This brings about another issue:
3. Reinitialization: After you made the discontinuous change to your model,
consistent initial conditions to continue simulating continuous time
behavior have to be found. This is the topic of what is discussed in the
quoted material in this missive.
The reinitialization can be implicit or explicit. For example, in case of a
bouncing ball, you would reinitialize its velocity when it hits the floor by
reverting its velocity upon impact, possibly accounting for a coefficient of
restitution (< 1).
In other cases, the consistent initial conditions that you need are more
difficult to obtain as they are implicit in the new model structure (this
happens for differential and algebraic equations [DAEs] of 'high index').
For example, if you close a switch (modeled as ideal) between two capacitors
connected in parallel to ground, there is an instantaneous redistribution of
their charge to make sure they have the same voltage drop.
The new charge distribution can actually be computed from the system of
equations such that conservation of charge is maintained. There is a fairly
extensive body of literature on this topic. A paper I like myself is
@ARTICLE{verghese81,
author = "George C. Verghese and Bernard C. L}vy and Thomas
Kailath",
title = "A Generalized State-Space for Singular Systems",
journal = "IEEE Transactions on Automatic Control",
volume = 26,
number = 4,
month = aug,
pages = "811--831",
year = 1981,
}
In the past, I have used a numerical routine by Andras Varga to compute the
consistent initial conditions by deriving a pseudo-Weierstrass form first.
That worked well.
Note that just a minimum norm projection will, in general, result in
physically incorrect values (i.e., you may generate charge or momentum or
any other extensive quantity)
4. Reduced dynamics: Once you have computed consistent initial conditions,
you have to make sure you keep generating consistent values. This is not
trivial in case your 'generalized state space' is of higher dimension than
your actual state space. In this case you live in a subspace and you want to
make sure that you do not leave this subspace during numerical integration
of the continuous time behavior.
There is quite some literature on this topic as well. An often cited
algorithm was presented by Pantelides in
@ARTICLE{pantelides88,
author = "Constantinos C. Pantelides",
title = "The Consistent Initialization of Differential-Algebraic
Systems",
journal = "SIAM Journal of Scientific and Statistical Computing",
volume = 9,
number = 2,
pages = "213--231",
month = mar,
year = 1988,
}
5. Sequences of mode changes: Now you can also have a number of
zero-crossing functions changing their sign as a result of the newly
computed state variable values. This causes a sequence of state changes,
with their corresponding complications. Let me just refer to a paper I wrote
http://moncs.cs.mcgill.ca/people/mosterman/papers/wintersim2003/p.pdf, which
gives a nice overview.
Anyway, I hope all this clarifies the complexity of simulating hybrid
dynamic systems a bit, or, at least, helps in structuring whatever
discussion is to follow.
Have a great weekend, y'all!
Pieter
--
Pieter J. Mosterman
Real-Time and Modeling & Simulation Technologies
The MathWorks, Inc.
http://www.mathworks.com
http://moncs.cs.mcgill.ca/people/mosterman

Yes.
The main problem with Laplace transforms is that they only work
for cases with constant coefficients (or sometimes when the
coefficients are simple functions). They are not exactly ideal
for the general case. Thus, e.g., if the original equation had
been
y'' + ( k^2 - A*\delta (x) -u(x) ) * y = 0
where u(x) is some continuous function, you would have a problem.
Moreover, I placed the discontinuity at x=0, assuming one wanted
the solution for both x>0 and x<0. One would either need the two-
sided Laplace transform or else to shift the point of discontinuity
as you did.
IIRC, the two-sided Laplace transform is discussed in A. Erdelyi's
little book, "Operational calculus and generalized functions"
(Holt, Rinehart and Winston, New York, 1962), as well as in
B. Van der Pol & H. Bremmer, "Operational calculus based on the
two-sided Laplace integral" (Cambridge U. Press, Cambridge, 1955).

Simulink can do this. The blocks with discontinuity have a 'zero
crossing detection' function that finds the location of the
discontinuity and then integrates on either side of it.
Don't know your problem but it may be a typical 'state event' type of
thing where you need a different structure before and after the event.
Simulink's enabled subsystems handles this nicely. Take at look
at the clutch demo.
dave y.

where
that
and
am
such
perhaps
If you really didn't want to buy any more software, you could use
runge-kutta
http://mathworld.wolfram.com/Runge-KuttaMethod.html
It may get messy, but you can definitely use discontinuous values here as
well. Note this is not a variable step length integrator.
tim

You can use RK for equations with discontinuous right-hand-side but
you should not. For discussion why, see the book "Solving Ordinary
Differential Equations I" by Hairer, Norsett and Wanner, Section
II.5 "Further Questions of practical computation - Discontinuous
equations". Some numerical experiments with RK and discontinuous
equations are presented, with conclusions: "... the results are
quite chaotic. In general, one looses one or two digits of precision
and the number of steps is very large".
A.L.

Polytechforum.com is a website by engineers for engineers. It is not affiliated with any of manufacturers or vendors discussed here.
All logos and trade names are the property of their respective owners.