Thread Subject: how to circumvent the removable singularities in ODE?

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!

Luna Moon <lunamoonmoon@gmail.com> wrote in message <70332e5b-9248-4dca-8048-6444fbdec6e4@f20g2000prn.googlegroups.com>...
> 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.
>

Take a look at stiff solver with mass matrix such as ode15s, ode23t.

Bruno
> 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!

In article <70332e5b-9248-4dca-8048-6444fbdec6e4@f20g2000prn.googlegroups.com>,
Luna Moon <lunamoonmoon@gmail.com> wrote:
>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!

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

On Nov 12, 3:58 pm, fakeu...@invalid.domain wrote:
> In article <70332e5b-9248-4dca-8048-6444fbdec...@f20g2000prn.googlegroups.com>,
> Luna Moon  <lunamoonm...@gmail.com> wrote:
>
>
>
> >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!
>
> 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

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