# Applying l'Hôpital's Rule to NDSolve with Piecewise Function

OK, so I've been working on this problem for over a week and it's driving me crazy.

I'm trying to solve a system of two differential equations in $x(t)$ and $y(t)$ using a "time-elimination method," such that $\tfrac{d x}{d t}/\tfrac{d y}{d t}=\tfrac{dx}{dy}$ reduces the system to a single differential equation in $x(y)$. The tricky part is that the initial condition needs to be defined at the steady state of the system $\{x_{ss},y_{ss}\}$, where $x'(t)=y'(t)=0$. The slope at this point is undefined since $dx/dy=0/0$ as $y\to y_{ss}$ by construction.

I've seen two common ways of dealing with this. The first involves finding a kind of approximation to the solution by slightly perturbing the correct initial condition so that $x'(y)$ is always well defined. For example, suppose the system is defined as:

y0 = 1; yss = 10;
xdot = 2*x[t]-y[t];
ydot = (.01 + .1*x[t])*(Log[yss/y[t]]);


We can find a solution using:

xdoty[y_] = (xdot/ydot) /. {x[t] -> x[y], y[t] -> y};
NDSolve[{x'[y] == xdoty[y], x[yss - .001] == xss}, x, {y, y0, yss - .001}];


which yields the following results for $x(t)$ (purple) and $y(t)$ (blue):

Notice that the solution for $y(t)$ "blows up" as it comes within .001 of $y_{ss}$.

The second, more accurate approach is to overcome the ill-defined slope at the steady state via an application of l'Hôpital's rule.

In Mathematica, I usually deal with this by passing a Piecewise function to NDSolve, e.g.:

NDSolve[{x'[y]==Piecewise[{{xdoty[y],y<yss},{xdotyss,y==yss}}],x[yss]==xss},x,{y,y0,yss}]


where xdotyss is $x'(y_{ss})$ (a single scalar value) calculated after the application of l'Hôpital's rule. This yields the following results for the system:

where both $x(t)$ and $y(t)$ are well-behaved throughout.

The problem I am having is that this approach does not always seem to work well. Particularly, for some differential systems I can get the first "perturbation" method to give me approximately correct looking results, but the second "l'Hôpital method" always gives me the following error:

NDSolve::ndsz: At y == yss, step size is effectively zero; singularity or stiff system suspected.

Given that the first approach causes no errors, and the second approach is using the same differential equations to derive the value of the slope at that single point, shouldn't this always be a viable approach? I can't think of what would be causing the system to suddenly become "stiff."

Any ideas?

-
Some of the variables don't seem to be defined. When I run the code for example, xss doesn't seem to be defined. Please try including a complete example. –  Searke Mar 15 '13 at 20:21
xss is just the value of x[t] that sets xdot=0 after substituting in yss for y[t]. I didn't include all the steps in the sample above, because it was not actually demonstrating the problem I was experiencing in another context. –  clr66 Mar 15 '13 at 23:28
Well, I guess I just needed one more day to work through the problem. Apparently using the Piecewise function as I outlined above does work generally. My issue with getting the ndsz "stiffness" error above apparently had something to do with precision, as changing all of my parameters to infinite precision forms (e.g. a=7/10 instead of a=0.7) allowed NDSolve to complete without any problems. Hooray! I hope someone else might find this useful. In any case, I'm very grateful to anyone who took a moment to look at the question and give it a thought. Cheers!