26
votes
Accepted
How can I avoid catastrophic cancellation?
Take
\begin{align}
1-\sqrt{ 1-x^2} &= (1-\sqrt{ 1-x^2})\frac{1+\sqrt{ 1-x^2}}{1+\sqrt{ 1-x^2}}\\
&= \frac{x^2}{1+\sqrt{ 1-x^2}}
\end{align}
So
\begin{align}
y = x\sqrt{\frac{1}{2+2\sqrt{1-x^2}}...
- 725
25
votes
Meaning of "-0.0" in Python?
Floating point numbers (according to the standard1 nearly all programming languages use) are stored with a certain number of bits in the mantissa, in the exponent, and with a sign bit. As such, ...
- 56.2k
24
votes
What language should I use when teaching an undergraduate course in computer programming?
In 2014, I would've said Python. In 2017, I wholeheartedly believe that the language to teach undergraduates is Julia.
Teaching is always about a tradeoff. On one hand, you want to choose something ...
- 12.4k
21
votes
How much more work is it to code math models in Python, compared to working with Matlab?
Going from MATLAB to Python does introduce quite a bit of syntax overhead. One way to quantify it is the nice QuantEcon cheatsheet which showcases how there's a lot of extra "stuff" going on ...
- 12.4k
17
votes
Accepted
CUDA & Python for numerical integration and solving differential equations
Julia's DifferentialEquations.jl is all GPU-compatible. If you make your arrays GPU-based arrays, then the solver recompiles to be all on the GPU (no data transfers). For example:
...
- 12.4k
15
votes
Tanh-sinh quadrature numerical integration method converging to wrong value
There are a bunch of pitfalls when it comes to tanh-sinh quadrature, one being that the integrand needs to be evaluated very closely to the interval boundaries, at distances less than machine ...
- 3,136
15
votes
Matrix multiplication accuracy Matlab vs Python
First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :)
Third, the ...
- 12k
15
votes
Time and memory required to diagonalize a 18000 by 18000 matrix using numpy in python
A 20000 by 20000 double-precision complex matrix requires
$20000 \times 20000 \times 8 \times 2=6.4 \mbox{gigabytes}$
of RAM. The LAPACK routines ZHEEV that will do the work for you will store the ...
- 18.9k
14
votes
Why does Matlab's integral outperform integrate.quad in Scipy?
The question has two very different subquestions. I will address the first one only.
Matlab's version runs on average 24 times faster than my python
equivalent!
The second one is subjective. I ...
- 141
14
votes
13
votes
How do I reliably generate random numbers in Python distributed across multiple nodes?
To the best of my knowledge, Numpy does not support independent streams. Indeed, getting independent streams from the Mersenne Twister (Pythons RNG) is notoriously difficult although it can be done.
...
- 2,513
13
votes
Accepted
Comparing Algorithmic complexity, ODE Solvers (Big O)
odeint from the SciPy library defaults to the lsoda integrator described here. However, any simple description of asymptotic ...
- 12.4k
13
votes
Accepted
How much more work is it to code math models in Python, compared to working with Matlab?
There are libraries that you can use in Python that will give you all (or at least nearly all) of the functionality of MATLAB. For example, scipy.integrate.solve_ivp() supports a number of methods ...
- 18.9k
12
votes
Runge-Kutta in the presence of an attractor
The problem you are encountering is likely not a consequence of your choice of algorithm, but in fact a consequence of the resulting dynamical system after applying time reversal. Per the definition ...
- 2,936
11
votes
Parallelizing a for-loop in Python
Without assuming something special on my_function choosing multiprocessing.Pool().map() is a good guess for parallelizing such ...
- 211
11
votes
Finding the first N roots of transcendental equation
This is a root-finding problem for an analytic function over a single dimension. One standard technique is to approximate $F(k)$ as a polynomial $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ using a ...
- 2,495
11
votes
Accepted
Arbitrary Precision Optimization Libraries?
Optim.jl from Julia will work with the number types that you give it, so if you make it use BigFloats then it'll do that. Local derivative based, derivative-free, global, and integrates with automatic ...
- 12.4k
11
votes
Accepted
Generate random smooth 2D closed curves
Since your figure is a closed loop, its parametric curves $x(t)$ and $y(t)$ must be periodic functions. This suggests one way to generate such figures, by constructing random smooth periodic functions ...
- 5,046
10
votes
Are there any "light-weight" FEM packages around?
I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you ...
- 2,104
10
votes
Accepted
Computing numeric derivative via FFT - SciPy
FFT returns a complex array that has the same dimensions as the input array. The output array is ordered as follows:
Element 0 contains the zero frequency component, F0.
The array element F1 ...
- 2,595
10
votes
How to solve a second order differential equation (diffusion) with boundary conditions using Python
I have found that I must keep the value of dt near dx or the results become unstable
This behavior you have noticed is known as the Courant–Friedrichs–Lewy (CFL) condition. Indeed, there are ...
- 1,124
9
votes
Is there a high quality nonlinear programming solver for Python?
We recently released (2018) the GEKKO Python package for nonlinear programming with solvers such as IPOPT, APOPT, BPOPT, MINOS, and SNOPT with active set and interior point methods. One of the issues ...
- 506
9
votes
Accepted
Correct eigenfunctions of Laplace operator by Finite Differences
You should specify the eigenvalues you want with which="SM", for example.
Check the following snippet. I also changed the solver, since your system is symmetric.
<...
- 8,582
9
votes
Accepted
Poor SVD reconstruction of singular matrix
Algorithms for the SVD, as more or less every classical linear algebra algorithm based on orthogonal transformations, are normwise backward stable, i.e., it should be guaranteed that $\frac{\|USV^* - ...
- 12k
9
votes
Accepted
Loop optimization with f2py, Cython and Numba
I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS ...
- 3,839
9
votes
Recommended language/environment for large scale semi-continuous biological models
You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq.
...
- 12.4k
9
votes
Accepted
How to choose a python parallelization library?
Dask schedules tasks across processes and across nodes, so it is appropriate for use on a single computer, supercomputer, or cloud. Dask also provides specialized data structures to aid in this.
...
- 4,021
9
votes
Computing numeric derivative via FFT - SciPy
Maxim Umansky’s answer describes the storage convention of the FFT frequency components in detail, but doesn’t necessarily explain why the original code didn’t work. There are three main problems in ...
- 191
9
votes
Accepted
Efficiently computing $e^{tX}$ for many different values of $t$
An anti-Hermitian matrix is diagonalizable, with orthogonal eigenvectors (ref). Hence you can write $X = PDP^{-1}$, where $D$ is a diagonal matrix. Therefore the exponential can be calculated as $e^X=...
- 254
8
votes
Accepted
Visually appealing ways to plot singular vector fields with matplotlib or other foss tools
[I took your sample program as a starting point and adapted Colormap Normalization from the matplotlib wiki.]
Almost everything of the picture just looks red.
Indeed. They problem is that there ...
- 465
Only top scored, non community-wiki answers of a minimum length are eligible