numerical-methods

Can someone explain why scipy.integrate.quad gives different results for equally long ranges while integrating sin(X)?

I am trying to numerically integrate an arbitrary (known when I code) function in my program using numerical integration methods. I am using Python 2.5.2 along with SciPy's numerical integration package. In order to get a feel for it, i decided to try integrating sin(x) and observed this behavior- >>> from math import pi >>> from scipy....

How to do numerical integration with quantum harmonic oscillator wavefunction?

How to do numerical integration (what numerical method, and what tricks to use) for one-dimensional integration over infinite range, where one or more functions in the integrand are 1d quantum harmonic oscillator wave functions. Among others I want to calculate matrix elements of some function in the harmonic oscillator basis: phin(x...

How can I use numpy.correlate to do autocorrelation?

Hi, 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 result, as it almost always gives a vector where the first number is not the largest, as it ought to be. So, this question is really tw...

Algorithm for multidimensional optimization / root-finding / something

I have five values, A, B, C, D and E. Given the constraint A + B + C + D + E = 1, and five functions F(A), F(B), F(C), F(D), F(E), I need to solve for A through E such that F(A) = F(B) = F(C) = F(D) = F(E). What's the best algorithm/approach to use for this? I don't care if I have to write it myself, I would just like to know where to...

How to make topographic map from sparse sampling data?

I need to make a topographic map of a terrain for which I have only fairly sparse samples of (x, y, altitude) data. Obviously I can't make a completely accurate map, but I would like one that is in some sense "smooth". I need to quantify "smoothness" (probably the reciprocal the average of the square of the surface curvature) and I wan...

implementing the derivative in C/C++

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

std::pow gives a wrong approximation for fractional exponents

Here is what I mean trying to do double x=1.1402 double pow=1/3; std::pow(x,pow) -1; result is 0 but I expect 0.4465 the equation is (1 + x) ^3= 1.1402, find x. ...

Function approximation

I have a function, P(x0, x1, ..., xn) that takes 100 integers as input and gives as output an integer. P is a slow function to evaluate(it can range from 30 seconds to a couple of minutes). I need to know which values of points will maximize the yielded value from P. What techniques can I use to accomplish this? I know generally peop...

Implementing tridiagonal matrix algorithm (TDMA) with NumPy

Hi. I'm implementing TDMA in Python using NumPy. The tridiagonal matrix is stored in three arrays: a = array([...]) b = array([...]) c = array([...]) I'd like to calculate alpha-coefficients efficiently. The algorithm is as follows: # n = size of the given matrix - 1 alpha = zeros(n) alpha[0] = b[0] / c[0] for i in range(n-1): a...

Using ARPACK with PARDISO

This is a question similar to a question here. I wonder is there already an open source implementation or example of ARPACK Eigensolver that works well with PARDISO solver and Intel mkl library? ...

Fastest numerical solution of a real cubic polynomial?

R question: Looking for the fastest way to NUMERICALLY solve a bunch of arbitrary cubics known to have real coeffs and three real roots. The polyroot function in R is reported to use Jenkins-Traub's algorithm 419 for complex polynomials, but for real polynomials the authors refer to their earlier work. What are the faster options for a r...

Problem with implementing successive approximation algorithm in C++ with Visual Studio 2008

I'm trying to implement an algorithm we did in numerical methods class. I have an identical procedure written in Maple and it works fine. I can't understand why it isn't working in C++. Any help or tips would be appreciated; thanks in advance. #include "stdafx.h" #include "math.h" #include <iostream> using namespace std; double absv...

positive semi-definite matrices and numerical stability ?

i'm trying to do factor analysis for a co-occurrence matrix(C) , which is computed from the term-document matrix(TD) as follows: C=TD*TD' In theory C should be positive semi-definite , but it isn't and the factor analysis algorithm can't work with it because of this.I can't change the algo because of speed reasons). I look it up and ...

Converting a repeating binary number to decimal (express as a series?)

Given a binary number that repeats, for example 0.(0011) or 0.0(101), how would one go about converting it to decimal? What I've been able to dig up so far is the simple method for converting a terminating binary number to decimal, as below: res(N+2) = res(N+1) / 2 + res(N) where res is the result after step N, and N is the current i...

Java library for simple numerical approximation

I need a Java library for numerical approximation. Nothing too complex. Just something that solves the following problem for example: epsilon = f(x) + y + z; find x such that epsilon gets close to zero. EDIT: f(x) is an arbitrarily complex function. It should accept parameters such as tolerance, etc. Aside from rolling it myself, j...

problem with arithmetic using logarthms to avoid numerical underflow

Hi! I have two lists of fractions; say A = [ 1/212, 5/212, 3/212, ... ] and B = [ 4/143, 7/143, 2/143, ... ]. If we define A' = a[0] * a[1] * a[2] * ... and B' = b[0] * b[1] * b[2] * ... I want to calculate the values of A' / B', My trouble is A are B are both quite long and each value is small so calculating the product causes...

Problem with arithmetic using logarithms to avoid numerical underflow (take 2)

I have two lists of fractions; say A = [ 1/212, 5/212, 3/212, ... ] and B = [ 4/143, 7/143, 2/143, ... ]. If we define A' = a[0] * a[1] * a[2] * ... and B' = b[0] * b[1] * b[2] * ... I want to calculate a normalised value of A' and B' ie specifically the values of A' / (A'+B') and B' / (A'+B') My trouble is A are B are both qui...

How do I obtain the machine epsilon in R?

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

C or Ada for engineering computations?

Hi,as an engineer I currently use C to write programs dealing with numerical methods. I like C as it's very fast. I don't want to move to C++ and I have been reading a bit about Ada which has some very good sides. I believe that much of the software in big industries have been or more correctly were written in Ada. I would like to know...

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 y=log(x,2) which is 215. x is always a power of 2. The problem seems to be really simple. All what is required is to find the position of...