Step Function: Numerical Precision Matters

[Goose]

Hi, Friends

How to numerically do with step function? I mean, this function is very sensitive to the precision. For example, f(x) is a step function. When x is greater than some very small positive number, \delta, the numerical answer is 1, which agrees with the expected value. But when x < \delta, the computer cannot tell it from 0, therefore the numerical result, f(x), equals to ZERO, which is not what I want.

Thanks for help

Goose

 

[CompuChip]

Hi Goose.

I don't really see what kind of step function you are using and what it is good for. Usually a step function can be defined as $\Theta_a(x) \equiv \Theta(x - a)$ where
$$\Theta(x) = \begin{cases} 0 & \text{ if } x < 0 \\ \frac12 & \text{ if } x = 0 \\ 1 & \text{ otherwise} \end{cases}$$
In numerical calculations this also works fine... where does delta come in?

 

[chiro]

If you want a continuous analytic solution you can try looking up articles on Fourier transforms dealing with the signum function. If you get a book dealing with various types of waveforms like sawtooth waves it will have a chance of containing that function (I can't remember it off the top of my head).

If you don't care about an analytic solution, just do what the poster described above and use some hard and fast boundary rules. The continuous analog exhibits the Gibbs phenomenon so you'll see weird effects around the boundaries but that's to be expected.

 

[Goose]

Thanks. My step function is the same as what CompuChip described. I know the exact analytical meaning of this function. But my question is how to define this function by Fortran or other computation language. In my research, I simply wrote subroutine using if statement as CoompuChip said. But the result is not desirable because computer doesn't know a very small number located in the neighbor of ZERO. For example, if $x = 10^{-50}$, the result is 0, not the expected value, 1. This difference affects my research a lot, and then my result always doesn't agree with my adviser's. Is there another way to construct step function?

Thanks again.

 

[CompuChip]

Ah, I see the problem now.
In that case I must leave answering your question to more qualified people. The only suggestions I can make are to a) find a package / routine / whatever it's called in Fortran which allows you to handle arbitrary floating point numbers and/or scientific notation; b) if you are only interested in values close to zero, scale everything by a factor 10^60 such that 10^{-50} is mapped to 10^{10}; c) check how your advisor did it.

But I hope people more experienced in numerical calculations will offer a good solution.

 

[Redbelly98]

I agree with CompuChip's suggestion to scale the numbers, if you can do so. Or choose a different system of units, so that you are not dealing with such small values.

That being said, here is a question for Goose: are you using double precision variables? They would handle numbers as small as 10^-307, or somewhere thereabouts, before treating the number as zero.

 

[Goose]

Good suggestions. Thank you all. I'll try it again.

 

[chiro]

Thanks. My step function is the same as what CompuChip described. I know the exact analytical meaning of this function. But my question is how to define this function by Fortran or other computation language. In my research, I simply wrote subroutine using if statement as CoompuChip said. But the result is not desirable because computer doesn't know a very small number located in the neighbor of ZERO. For example, if $x = 10^{-50}$, the result is 0, not the expected value, 1. This difference affects my research a lot, and then my result always doesn't agree with my adviser's. Is there another way to construct step function?

Thanks again.

If your using the explicit branching definition then you could probably do one of two things:

a) Define an absolute minimum delta value and always compare against that. Note that anything between the delta and zero will naturally using this definition be defined as zero and not as the number with regards to function behaviour.

b) Compare the floats memory specific value to a hardcoded value of 0. This is a few extra lines because of the intrinsics of how floating point numbers are represented with the mantissa and with things such as NaN and so forth. One hint I can say from experience is that -0 and +0 are both zero so you should always add +0 to the number before checking the footprint to avoid that ambiguity.

You asked if there was any other way to do this. The best way is to do what you've done and use an explicit branch condition. Another way is to use the analytical version. It however is prone to so much error that you would need to know a bit about numeric analysis to make the output stable at certain values.

Anyway you've probably got it by now so i'll leave it there