11
votes
What algorithm(s) do numpy and scipy use to calculate matrix inverses?
Documentation to numpy.linalg.inv and scipy.linalg.inv does not mention the algorithm used.
Judging from the source, ...
- 325
11
votes
10th-order Runge-Kutta Method
Julia's DifferentialEquations.jl package (linked here via the subpackage OrdinaryDiffEq.jl which contains the source code of many of the methods) includes an implementation of a 10th order explicit RK ...
- 2,936
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 are scientific computing workflows faring on Apple's M1 hardware
The way I have been measuring whether the eco-system is ready is how things are going with the transition for the homebrew package manager. They have been carefully documenting the progress of ...
- 556
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
Is LAPACK behind the cutting edge of dense linear algebra?
When one says an algorithm is of order $O(n)$, that may mean that the complexity is given by: $c + b*n$. With every new element you add you increase in runtime (effectively). What mathematically ...
- 3,005
9
votes
Accepted
Number of function calls and jacobian calls in scipy.root
At the point where you print out the Jacobian, adding
traceback.print_stack()
reveals that the first evaluation comes from ...
- 716
9
votes
10th-order Runge-Kutta Method
You seem to be misunderstanding. Since you're referring to an explicit Runge-Kutta method, applying the method does not require that you solve any equations.
In Ernst Hairer's article, he is deriving ...
- 16.7k
8
votes
Accepted
Why does LSODA fail to integrate the logistic function?
When you use $r=5$, the initial condition is
$$ x(-10) \approx e^{-50} \approx 1.9\times 10^{-22}. $$
This is much smaller than the machine epsilon, $2\times 10^{-16}$, and it is very likely that ...
- 11.5k
8
votes
Is LAPACK behind the cutting edge of dense linear algebra?
LAPACK has been on the cutting edge for just about three decades, and probably still is for its niche. However, given given recent developments in libraries for the simpler BLAS-type matrix operations ...
- 181
8
votes
Accepted
Python scipy eigh(Arpack) giving wrong eigenvalues for generalized eigenvalue problem
The matrix B (M in the documentation) needs to positive definite according to the documentation: "If sigma is None, M is positive definite", this is in addition to the first requirement &...
- 407
8
votes
How are scientific computing workflows faring on Apple's M1 hardware
From this, it looks like there is no functional native Fortran compiler yet. If that is really the case, things look bleak. Almost anything that uses linear algebra includes some Fortran code (Lapack),...
- 12k
7
votes
Accepted
Matrix Balancing Algorithm
Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit.
After checking the problematic cases reported in David S. Watkins. A case where balancing is ...
- 393
7
votes
Accepted
Applying the result of Cuthill-McKee in SciPy
The reverse Cuthill-McKee algorithm produces a reordering that applies to both the rows and columns. This is because it works by considering matrices as graphs of (undirected) connected nodes. ...
- 188
7
votes
Accepted
Integration of the Fermi distribution using Python
One of your problems is the system of units that you are using. Just changing the units improves the results
...
- 8,582
7
votes
Accepted
Forward and backward integration -- cause of errors
Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale ...
- 12.4k
7
votes
Is there an efficient way to form this block matrix with numpy or scipy?
The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $A^i B$, with $i=0\,\dots,N$ you do not have to compute $A^i$ at each step, but you ...
- 3,839
7
votes
Accepted
Algorithm to factorize matrix whose many rows are already of upper triangular form?
I believe you can accomplish what you want efficiently using the recursive LU algorithm. In brief, recursive LU on a $M \times N$ matrix $A$ proceeds by partitioning the matrix into 4 blocks:
\begin{...
- 1,068
7
votes
Accepted
Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix
I think you are overestimating the overhead of computing L. There are zero extra operations needed; the only additional cost is writing to RAM some numbers that you ...
- 12k
7
votes
Accepted
Why is my curve_fit not producing the covariance matrix and the correct values for the unknown variables?
The problem seems to be one of scaling. When I added the jacobian of the function an overflow warning appeared. Thus, I divided the data by their maximum values and it worked. Following is the code.
<...
- 8,582
6
votes
Accepted
Fast c++ library to solve very big sparse systems
I second the idea of using Eigen, which is pretty efficient, but also very simple to include.
If you need a lot more performance, you could try to use PETSc or Trilinos. They are very powerful ...
- 1,157
6
votes
Accepted
Correct use of scipy's sparse.linalg.spilu
If $\mathbf L$ and $\mathbf U$ give an approximate factorization of $\mathbf A$, you wouldn't want to use $\mathbf P = \mathbf L\cdot \mathbf U$ as a preconditioner (that's approximately $\mathbf A$), ...
- 5,046
6
votes
Accepted
Solving for a set of coupled ODEs to get correct variable values
The function $q(e)$ satisfies a first order linear ODE
$$ \frac{\mathrm{d}q}{\mathrm{d}e} = \frac{111 e^4+876 e^2+288}{(e^2-1)
(121 e^2+304)} q(e), $$
which can be solved very easily by using an ...
- 11.5k
6
votes
Accepted
Ways to solve $Ax=b$ for a sparse (banded) $A$ with updates
If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1.
More generally, you can talk about matrices of bandwidth $k$ where $k$ is any ...
- 10.4k
6
votes
Numerical integration in Python with unknown constant
What you want seems inherently impossible, and that’s not due to restrictions of Python.
The only way we can arrive at a situation where we only need to apply a single quadrature is to get ...
- 2,127
6
votes
Accepted
Get the roots of a Hermite interpolating polynomial
The interpolated polynomial does not have roots. Considering that the behavior outside the interpolation region holds is termed extrapolation.
You can explicitly use the polynomial, given by (as I ...
- 8,582
6
votes
Numerical derivative in python
The Savitzky-Golay filter uses a constant delta (the spacing of the samples,) and the default value of the delta in the filter implementation is 1, according to https://docs.scipy.org/doc/scipy-0.16.1/...
- 81
6
votes
solve_ivp from scipy does not integrate the whole range of tspan
You can examine the sol object to see why the integration failed. It provides the message 'Required step size is less than spacing between numbers.' This usually ...
- 1,124
6
votes
Accepted
Calculate determinant of unitary matrices based on SVD implementation
If you are prepared to go digging around in the fortran code:
The SVD algorithm consists of a few parts:
Bidiagonalization (usually using Householder reflectors)
Use QR shifts to reduce the ...
- 1,785
Only top scored, non community-wiki answers of a minimum length are eligible