34
First, if your undergraduates are like ours and had no prior introduction to computers, expect to spend some time teaching them how to use basic stuff like using a proper editor (i.e., not MS Word), the command line, etc.
I think the answer somewhat depends on where you set the focus of your course (or what you are required to teach). For example: How ...
32
Ease of learning
Python and Fortran are both relatively easy-to-learn languages. It's probably easier to find good Python learning materials than good Fortran learning materials because Python is used more widely, and Fortran is currently considered a "specialty" language for numerical computing.
I believe the transition from Python to Fortran ...
31
This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$
Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that ...
25
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}}}
\end{align}
24
I am not super familiar with f2py internals, but I am very familiar with wrapping Fortran. F2py just automates some or all of the things below.
You first need to export to C using the iso_c_binding module, as described for example here:
http://fortran90.org/src/best-practices.html#interfacing-with-c
Disclaimer: I am the main author of the fortran90.org ...
22
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 that is simple enough that it is easy to grasp. But secondly, you want to teach something that has staying power, i.e. something that can grow with you. The ...
22
Joblib does what you want.
The basic usage pattern is:
from joblib import Parallel, delayed
def myfun(arg):
do_stuff
return result
results = Parallel(n_jobs=-1, verbose=verbosity_level, backend="threading")(
map(delayed(myfun), arg_instances))
where arg_instances is list of values for which myfun is computed in parallel.
The main ...
20
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 when trying to write simple linear algebra commands in Python. The verbose NumPy syntax is really just a symptom of how it was not developed as a technical ...
18
There are two issues that you are likely to be encountering.
Ill-conditioning
First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian:
$$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$
If you ...
17
Here is the Numba solution. On my machine the Numba version is >1000x faster than the python version without the decorator (for a 200x200 matrix, 'k' and 200-length vector 'a'). You can also use the @autojit decorator which adds about 10 microseconds per call so that the same code will work with multiple types.
from numba import jit, autojit
@jit('f8[:]...
16
For the first part of my question, I found this very useful comparison for performance of different linear interpolation methods using python libraries:
http://nbviewer.ipython.org/github/pierre-haessig/stodynprog/blob/master/stodynprog/linear_interp_benchmark.ipynb
Below is list of methods collected so far.
Standart interpolation, structured grid:
http:/...
15
A difficulty with any of these types of questions is that the answer is highly community-dependent.
To answer some of your questions in haphazard order:
MATLAB is used a lot both in academia and in industry. One of the reasons it's used quite a bit in industry is because it is taught in academia. I know for a fact that MATLAB is used at Lincoln Laboratory ...
15
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 would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...
14
What you're looking for is Numba, which can auto parallelize a for loop.
From their documentation
from numba import jit, prange
@jit
def parallel_sum(A):
sum = 0.0
for i in prange(A.shape[0]):
sum += A[i]
return sum
14
Here is R1, as computed in MATLAB:
1.0e+07 *
-7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830
-1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472
3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415
-5....
14
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 real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...
13
I will address only the comparison of C to C++. While it is true that anything written in C can be ported to C++ with a few syntactic touch-ups, the communities have different values.
The C library community, more than almost any other, values binary stability.
Binary stability is critical for low-level libraries to avoid inflicting constant pain on those ...
13
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.
Consider using the RandomGen package. It is fully compatible with Numpy, and provides you with the PCG64 generator, supporting up to $2^{63}$ independent ...
13
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 of an attractor, all points in some neighborhood of the attractor will converge to the attractor under the flow of the dynamical system as $t\to\infty$. However, ...
13
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 for ODE integration that are comparable to what you can get with the various odexxx() functions in MATLAB.
So no, you wouldn't have to write your own ODE ...
12
The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library.
In exact arithmetic, the Lanczos algorithm ...
11
According to the docs, there is no in-place permutation method in numpy, something like ndarray.sort.
So your options are (assuming that M is a $N\times N$ matrix and p the permutation vector)
implementing your own algorithm in C as an extension module (but in-place algorithms are hard, at least for me!)
$N$ memory overhead
for i in range(N):
M[:,i] = ...
11
CVXOPT only solves (smooth and nonsmooth) convex problems, giving access to several third party convex solvers with guaranteed state of the art worst case complexity. You may pose linear, convex quadratic, linear semidefinite, and many other convex types of constraints.
OpenOpt solves general (smooth and nonsmooth) nonlinear programs, including problems ...
11
As Misha and Geoff Oxberry pointed out, Mathematica really has a different focus (just because you can pound in a nail with a screwdriver doesn't mean you should). So I take your question as being "If I know Matlab, why should I learn Python?"
[Edit: and so, apparently, did you.]
For all intents and purposes, Matlab is the English of scientific computing -- ...
11
odeint from the SciPy library defaults to the lsoda integrator described here. However, any simple description of asymptotic computation time is impossible. The reason is many fold.
First, let me describe the algorithm. A common multistep algorithm for non-stiff equations are the Adams-Moulton methods. While these are implicit, the Adams-Bashforth methods ...
11
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 Chebyshev approximation, and to compute the roots of the polynomial semi-analytically, e.g. by setting up the companion matrix and computing its eigenvalues. (...
10
I have a suspicion that your Python expression for the right-hand side is doing integer division, not floating-point division. As a result, ((l1-l2)**2)/(d**2) is being evaluated as zero, and the term inside the square root is one.
In fact, you forgot the reciprocal in your right-hand expression as well, but that's not the first problem...
In MATLAB:
>&...
10
I was pondering this a few days ago (also in Python). Personally I don't think that object oriented programming is always a good fit for numerical programming. You can get distracted with designing the classes rather than just solving the equations. I prefer to stay with simple functions, and with numpy you can have your equations vectorised so the number of ...
10
Python is a very slow, high level language. For fast number crunching you'll have to write the main compute kernels in low level languages like C/C++ which means that now you have to learn not one but at least two languages. You'll also have to deal with additional headache associated with debugging/installation/maintenance etc. Most people use Python as a ...
10
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 precision, e.g., 1.0 - 1.0e-20 in the original example. When this point is evaluated, it rounds to 1.0 at which f has a singularity, and anything can happen. You're ...
Only top voted, non community-wiki answers of a minimum length are eligible
