how to circumvent the removable singularities in ODE?

[Luna Moon]

Hi all,

I am implementing an ODE in Matlab.

y'(x)=f(x),

for x in [a, b].

However, the f(x) is very complicated, there is a numerator and a
denominator,

both can be infinity or zero, because of the involvement of "log" and
"exp" in f(x).

However, when I use Maple and calculate the value of f(x) at
singularity points,

I found they are removable singularity points.

So if Matlab numerical evalution of f(x) is as smart as the symbolic
calculation of Maple, then it should be worry-free.

But it's not because of the numerical explosion near those singular
points.

Of course I could set a small threshold, and whenever x enters into
that small neighborhood, I switch to the analytical results obtained
from Maple limiting operation. However, choosing that threshold leads
to discontinuity, isn't it?

So what's the best way to handle such situation? I would like to hear
your thoughts. Thanks a lot!

[fake...@invalid.domain]

In article <70332e5b-9248-4dca...@f20g2000prn.googlegroups.com>,

Consider sin(x)/x with its removable singularity at zero. I'm not
sure about matlab, but scilab isn't smart enough on its own; there
is a function sinc(x) with the singularity removed.

I seem to remember that if then statements inside functions caused
problems when the inputs were vectors.

The scilab code below makes f2 to be f with the singularity removed.
the ode code fails with f but suceeds with f2.

function y = f(t,x)
y = sin(x)/x;
endfunction

function y = f2(t,x)
top = sin(x);
bottom = x;
top(find(x==0))=1;
bottom(find(x==0))=1;
y = top/bottom
endfunction

t=0:0.1:1;
w0=0;t0=0;
w =ode(w0,t0,t,f2)

matlab with the symbolic toolbox can call maple (at least in older
versions. Might now be mumath?

--
Steven Bellenot http://www.math.fsu.edu/~bellenot
Professor and Associate Chair phone: (850) 644-7405
Department of Mathematics office: 223 Love
Florida State University email: bellenot at math.fsu.edu

[RRogers]

On Nov 12, 3:58 pm, fakeu...@invalid.domain wrote:
> In article <70332e5b-9248-4dca-8048-6444fbdec...@f20g2000prn.googlegroups.com>,

As I recall I found a formulation/process in Maple that prompted the
numerical routines to invoke L'Hopital's rule. I am not sure but I
can dig up the code if you like. I think it was just a minor trick;
and probably gross.
Ray

[Helmut Jarausch]

Here some thoughts.
If you have (quite) a recent version of Matlab you probably can
use exception handling
http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_prog/f3-38012.html
(but I have no experience with it)

If you choose a threshold for both the size of the numerator and the size
of the denominator you can estimate the size of the tiny jump.
In the worst case the integrator notices it and reduces the step size
(unnecessarily). One option in that case is to stop the integration,
and restart it at a slightly greater time. The size of the jump is then
multiplied with the size of the time increment which should be really tiny.

--
Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany