# Questions tagged [numerical-methods]

Algorithms which solve mathematical problems by means of numerical approximation (as opposed to symbolic computation).

2,099
questions

125votes

13answers

226kviews

### How can I use numpy.correlate to do autocorrelation?

I need to do auto-correlation of a set of numbers, which as I understand it is just the correlation of the set with itself.
I've tried it using numpy's correlate function, but I don't believe the ...

118votes

5answers

239kviews

### Dealing with float precision in Javascript [duplicate]

I have a large amount of numeric values y in javascript. I want to group them by rounding them down to the nearest multiple of x and convert the result to a string.
How do I get around the annoying ...

90votes

1answer

8kviews

### Is there special significance to 16331239353195370.0?

Using import numpy as np I've noticed that
np.tan(np.pi/2)
gives the number in the title and not np.inf
16331239353195370.0
I'm curious about this number. Is it related to some system machine ...

65votes

11answers

142kviews

### Check if a varchar is a number (TSQL)

is there an easy way to figure out if a varchar is a number?
Examples:
abc123 --> no number
123 --> yes, its a number
Thanks :)

58votes

3answers

5kviews

### BigDecimal multiply by zero

I am performing a simple multiplication with BigDecimal and I have found some strange behaviour when multiplying by zero (multiplying by zero is correct in this use-case).
Basic maths tells me that ...

45votes

3answers

24kviews

### How do I obtain the machine epsilon in R?

Is there a constant that stores the machine epsilon in R?

37votes

5answers

63kviews

### Signal processing library in Java? [closed]

I'd like to compute power spectral density of time series; do some bandpass, lowpass, and highpass filtering; maybe some other basic stuff.
Is there a nice open-source Java library to do this?
I've ...

32votes

15answers

35kviews

### Sampling a random subset from an array

What is a clean way of taking a random sample, without replacement from an array in javascript? So suppose there is an array
x = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
and I want to randomly sample 5 ...

25votes

6answers

9kviews

### What's the numerically best way to calculate the average

what's the best way to calculate the average? With this question I want to know which algorithm for calculating the average is the best in a numerical sense. It should have the least rounding errors, ...

24votes

8answers

65kviews

### Implementing the derivative in C/C++

How is the derivative of a f(x) typically calculated programmatically to ensure maximum accuracy?
I am implementing the Newton-Raphson method, and it requires taking of the derivative of a function.

24votes

1answer

14kviews

### Fast algorithm to calculate Pi in parallel

I am starting to learn CUDA and I think calculating long digits of pi would be a nice, introductory project.
I have already implemented the simple Monte Carlo method which is easily parallelize-able. ...

22votes

12answers

5kviews

### As Java is to Scala, C++ is to ...?

Although C++0x is quite an improvement to C++ (type inference, anonymous functions, and so on), I have to say that Scala seems even better. The thing is that Scala only runs on the JVM, although it ...

21votes

2answers

25kviews

### __builtin_prefetch, How much does it read?

I'm trying to optimize some C++ (RK4) by using
__builtin_prefetch
I can't figure out how to prefetch a whole structure.
I don't understand how much of the const void *addr is read. I want to have ...

20votes

2answers

658views

### What algorithm is R using to calculate mean?

I am curious to know what algorithm R's mean function uses. Is there some reference to the numerical properties of this algorithm?
I found the following C code in summary.c:do_summary():
case ...

20votes

6answers

9kviews

### Integer cube root

I'm looking for fast code for 64-bit (unsigned) cube roots. (I'm using C and compiling with gcc, but I imagine most of the work required will be language- and compiler-agnostic.) I will denote by ...

18votes

8answers

31kviews

### How to find a binary logarithm very fast? (O(1) at best)

Is there any very fast method to find a binary logarithm of an integer number? For example, given a number
x=52656145834278593348959013841835216159447547700274555627155488768 such algorithm must find ...

18votes

1answer

12kviews

### sample random point in triangle [closed]

Suppose you have an arbitrary triangle with vertices A, B, and C. This paper (section 4.2) says that you can generate a random point, P, uniformly from within triangle ABC by the following convex ...

17votes

5answers

96kviews

### Making C code plot a graph automatically

I have written a program which writes a list of data to a '.dat' file with the intention of then plotting it separately using gnuplot. Is there a way of making my code plot it automatically? My ...

17votes

2answers

18kviews

### How can I plot the derivative of a graph in gnuplot?

I have a set of measurements of a variable over time. I have these measurements in a file called "results" with this format:
# time sample
0 5
12 43
234 342
etc...
I can easily plot ...

16votes

3answers

8kviews

### Why is Matlab's inv slow and inaccurate?

I read at a few places (in the doc and in this blog post : http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/ ) that the use of inv in Matlab is not recommended because it is slow and ...

16votes

2answers

22kviews

### What's the best way to calculate a numerical derivative in MATLAB?

(Note: This is intended to be a community Wiki.)
Suppose I have a set of points xi = {x0,x1,x2,...xn} and corresponding function values fi = f(xi) = {f0,f1,f2,...,fn}, where f(x) is, in general, an ...

16votes

2answers

14kviews

### How can I increase the number of subdivisions for functions in `scipy.integrate.dblquad`?

I'm using scipy.integrate.dblquad, and I get this error:
UserWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement ...
I want to increase ...

16votes

1answer

563views

### Emulating fixed precision in python

For a university course in numerical analysis we are transitioning from Maple to a combination of Numpy and Sympy for various illustrations of the course material. This is because the students already ...

16votes

0answers

1kviews

### How to call UMFPACK as MATLAB does?

The problem
I wish to solve a general system of linear equations A*x=b. The m-by-m matrix is sparse, real, square, somewhat poorly-conditioned, non-symmetric, but it is singular (rank(A)==m-1) ...

15votes

4answers

21kviews

### Why is the output of inv() and pinv() not equal in Matlab and Octave?

I have noticed that if A is a NxN matrix and it has the inverse matrix. But what the inv() and pinv() function output is different.
- My environment is Win7x64 SP1, Matlab R2012a, Cygwin Octave 3.6.4,...

15votes

9answers

9kviews

### Generating digits of square root of 2

I want to generate the digits of the square root of two to 3 million digits.
I am aware of Newton-Raphson but I don't have much clue how to implement it in C or C++ due to lack of biginteger support. ...

15votes

5answers

6kviews

### Where do I find the machine epsilon in C#?

The machine epsilon is canonically defined as the smallest number which added to one, gives a result different from one.
There is a Double.Epsilon but the name is very misleading: it is the smallest (...

14votes

4answers

54kviews

### Numerical ODE solving in Python

How do I numerically solve an ODE in Python?
Consider
\ddot{u}(\phi) = -u + \sqrt{u}
with the following conditions
u(0) = 1.49907
and
\dot{u}(0) = 0
with the constraint
0 <= \phi <= 7\...

14votes

2answers

15kviews

### is there a c++ library for ordinary differential equation (ODE) solvers?

More specifically, i'm interested in 8th order Dormand-Prince embedded method, it's based on Runge-Kutta, and stiff equations.
I use Numerical Recipes 3 but i often have trouble compiling their ...

14votes

5answers

6kviews

### efficiently determining if a polynomial has a root in the interval [0,T]

I have polynomials of nontrivial degree (4+) and need to robustly and efficiently determine whether or not they have a root in the interval [0,T]. The precise location or number of roots don't concern ...

14votes

2answers

5kviews

### N body simulation in C++

I am trying to implement an OpenMP version of the 2-dimensional n-body simulation.
But there is a problem: I assume each particle's initial velocity and acceleration are zero. When the particles ...

14votes

4answers

1kviews

### Which numerical library to use for porting from Matlab to C++? [closed]

I am currently prototyping some algorithms in Matlab that rely on matrix, DSP, statistics and image analysis functionality.
Some examples of what I may need:
eigenvectors
convolution in 2D and 3D
...

14votes

2answers

272views

### Which floating-point comparison is more accurate, and why?

I'm experimenting with different implementations of the Newton method for calculating square roots. One important decision is when to terminate the algorithm.
Obviously it won't do to use the ...

14votes

2answers

1kviews

### Solving PDE with implicit euler in python - incorrect output

I will try and explain exactly what's going on and my issue.
This is a bit mathy and SO doesn't support latex, so sadly I had to resort to images. I hope that's okay.
I don't know why it's inverted, ...

14votes

1answer

5kviews

### Calculate eigenvalues/eigenvectors of hundreds of small matrices using CUDA

I have a question on the eigen-decomposition of hundreds of small matrices using CUDA.
I need to calculate the eigenvalues and eigenvectors of hundreds (e.g. 500) of small (64-by-64) real symmetric ...

13votes

4answers

9kviews

### argsort for a multidimensional ndarray

I'm trying to get the indices to sort a multidimensional array by the last axis, e.g.
>>> a = np.array([[3,1,2],[8,9,2]])
And I'd like indices i such that,
>>> a[i]
array([[1, 2, ...

13votes

7answers

3kviews

### Automatic differentiation library in Scheme / Common Lisp / Clojure

I've heard that one of McCarthy's original motivations for inventing Lisp was to write a system for automatic differentiation. Despite this, my Google searches haven't yielded any libraries/macros for ...

13votes

6answers

2kviews

### Numerical library for Scala

I'm looking for a library to do numerical computing in Scala (or Java, although something that can use scala functions would be way nicer!) with at least the following capabilities:
L-BFGS
Minimizers ...

13votes

2answers

3kviews

### Bezier Cubic Curves: moving with uniform acceleration

Let's say I have a Bezier curve B(u), if I increment u parameter at a constant rate I don't obtain a costant speed movement along the curve, because the relation between u parameter and the point ...

13votes

0answers

845views

### Is boost::math::sinc_pi unnecessarily complicated?

This is not a question about template hacks or dealing with compiler quirks. I understand why the Boost libraries are the way they are. This is about the actual algorithm used for the sinc_pi function ...

12votes

1answer

21kviews

### numpy gradient function and numerical derivatives

The array the numpy.gradient function returns depends on the number of data-points/spacing of the data-points. Is this expected behaviour? For example:
y = lambda x: x
x1 = np.arange(0,10,1)
x2 = np....

12votes

1answer

2kviews

### Binary GCD Algorithm vs. Euclid's Algorithm on modern computers

http://en.wikipedia.org/wiki/Binary_GCD_algorithm
This Wikipedia entry has a very dissatisfying implication: the Binary GCD algorithm was at one time as much as 60% more efficient than the standard ...

12votes

2answers

824views

### Why are structs so much faster than classes for this specific case?

I have three cases to test the relative performance of classes, classes with inheritence and structs. These are to be used for tight loops so performance counts. Dot products are used as part of many ...

11votes

3answers

7kviews

### how can I evaluate the derivative of a spline function in R?

R can generate a spline function using splinefun() in the splines library. However, I need to evaluate this function at its first and second derivatives. Is there a way to do this?
for example
...

11votes

6answers

3kviews

### How do I generate points that match a histogram?

I am working on a simulation system. I will soon have experimental data (histograms) for the real-world distribution of values for several simulation inputs.
When the simulation runs, I would ...

11votes

2answers

17kviews

### Solving equation using bisection method

Is there a bisection method I can find online, specifically for python?
For example, given these equations how can I solve them using the bisection method?
x^3 = 9
3 * x^3 + x^2 = x + 5
cos^2x +...

11votes

3answers

5kviews

### Robust polygon normal calculation

is there a good robust algorithm to calculate normal vector of a convex polygon (in 3D, of course)? For triangles, it is easy: one takes two of the triangle's edges and calculates the cross product:
...

11votes

1answer

314views

### Recursive arc-length reparameterization of an arbitrary curve

I have a 3D parametric curve defined as P(t) = [x(t), y(t), z(t)].
I'm looking for a function to reparametrize this curve in terms of arc-length. I'm using OpenSCAD, which is a declarative language ...

10votes

2answers

4kviews

### Numerically stable method for solving quadratic equations

Using floating point, It is known that the quadratic formula does not work well for b^2>>4ac, because it will produce a loss of significance, as it is explained here.
I am asked to find a better way ...

10votes

6answers

6kviews

### Generating Full Period/Full Cycle Random Numbers or Permutations Similar to LCG but without odd/even

I wish to generate psuedo-random numbers/permutations that 'occupy' a full period or full cycle within a range. Usually an 'Linear Congruential Generator' (LCG) can be used to generate such sequences, ...