Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
numerical error?

Subject: numerical error?

From: Alex James

Date: 19 Mar, 2010 22:40:09

Message: 1 of 4

I've run into an obnoxious problem with the below function behavior around zero. Not only does it have problems AT zero due to the divide by zero, it also has problems approaching zero, even for smallish values of nn (eg 11 or 101), which makes patching the divide by zero difficult as I'm not sure how many points surrounding zero to remove.

% fChandra.m - Chandrasekhar function as described in Smith PoP 2005
function fChandra=fChandra(xx) % the Chandrasekhar function used in Dreicer, Phys Rev 1959, eqn 20

    fChandra=(erf(xx)-xx.*gradient(erf(xx),xx))./(xx.*xx);
    fChandra( find(xx(:)==0) )=0.0; % fix divide by zeros

    % the function has numerical problems at zero, which the above fix doesn't entirely solve...
    %{
    nn=1e4+1;
    xx=linspace(-5,5,nn);find(xx==0)
    figure, hold all
    plot(xx,fChandra(xx))
    plot(xx,gradient(fChandra(xx),xx))
    plot(xx,erf(xx))
    plot(xx,-xx.*gradient(erf(xx),xx))
    plot(xx,1./(xx.^2))
    legend('chand','grad(chand)','erf','-xx*grad(erf,xx)','1/xx^2')
    ylim([-1.5,1.5])
    %}

Subject: numerical error?

From: Rune Allnor

Date: 19 Mar, 2010 23:37:45

Message: 2 of 4

On 19 Mar, 23:40, "Alex James" <alex.n.ja...@gmail.com> wrote:
> I've run into an obnoxious problem with the below function behavior around zero. Not only does it have problems AT zero due to the divide by zero

>> xx=0;
>> fChandra=(erf(xx)-xx.*gradient(erf(xx),xx))

fChandra =

     0

Since you end up with a 0/0 type expression, you ought
to be able to handle this using L'Hopital's rule. The
problems in the vicinity of x=0 are likely caused by
the numerical gradient. If you can find an analytical
expression for the gradient (which you will need anyway,
to employ L'Hopital's rule) the problems might go away.

Rune

Subject: numerical error?

From: Jan Simon

Date: 19 Mar, 2010 23:47:05

Message: 3 of 4

Dear Alex!

> fChandra=(erf(xx)-xx.*gradient(erf(xx),xx))./(xx.*xx);
> fChandra( find(xx(:)==0) )=0.0; % fix divide by zeros

Instead of gradient(erf(xx)) you can look at:
  http://mathworld.wolfram.com/Erf.html
to get the derivative:
  (2 / pi) * exp(-xx .* xx)

BTW: You do not need a FIND here:
  fChandra( find(xx(:)==0) )=0.0
==> faster
  fChandra((xx(:)==0) = 0.0
Good luck, Jan

Subject: numerical error?

From: Roger Stafford

Date: 20 Mar, 2010 00:40:09

Message: 4 of 4

"Alex James" <alex.n.james@gmail.com> wrote in message <ho0ug9$6m7$1@fred.mathworks.com>...
> I've run into an obnoxious problem with the below function behavior around zero. Not only does it have problems AT zero due to the divide by zero, it also has problems approaching zero, even for smallish values of nn (eg 11 or 101), which makes patching the divide by zero difficult as I'm not sure how many points surrounding zero to remove.
>
> % fChandra.m - Chandrasekhar function as described in Smith PoP 2005
> function fChandra=fChandra(xx) % the Chandrasekhar function used in Dreicer, Phys Rev 1959, eqn 20
>
> fChandra=(erf(xx)-xx.*gradient(erf(xx),xx))./(xx.*xx);
> fChandra( find(xx(:)==0) )=0.0; % fix divide by zeros
>
> % the function has numerical problems at zero, which the above fix doesn't entirely solve...

  As Rune and Jan have indicated, the use of matlab's gradient function is totally inappropriate here when a simple analytic expression for the derivative exists:

 d erf(x) / dx = 2/sqrt(pi)*exp(-x^2)

The gradient function is a numerical approximation which uses derivatives of polynomial approximations of functions based on the intervals at which the function is evaluated. This can be very inaccurate for a problem of your kind.

However, even when you have made that replacement, the form of your function is subject to errors in the vicinity of x = 0 because both the numerator and denominator approach zero as x approaches zero, and worse, the numerator is the difference between two terms which are much larger than their difference, which is an anathema to accurate computing As Rune suggests, in the immediate neighborhood of x = 0, the function values should be replaces by the first few terms of a taylor expansion, which if my computations are correct are:

 2/sqrt(pi)*(2/3*x - 2/5*x^3 + 1/7*x^5 - ...)

  Doing something of this kind is the only way you are going to avoid trouble in the neighborhood of x = 0.

Roger Stafford

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us