Hot answers tagged

13 votes
Accepted

Forcing an ODE solver to preserve the norm

The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are ...
user avatar
11 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 ...
user avatar
9 votes
Accepted

Algorithm for high quality 1/f noise?

You can generate a noise sequence with whatever noise spectrum you want (including $1/f$, also known as pink noise) by generating the noise coefficients in spectral space. The magnitudes of the ...
user avatar
9 votes

Methods for solving $x'=Ax+b$ for small, sparse, singular $A$

Any general-purpose ODE solver should be able to handle this linear coupled system of ODE very easily, for example: scipy.integrate.ode CVODE from the Sundials solver suite; it appears to have Python ...
user avatar
  • 11.4k
9 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),...
user avatar
8 votes
Accepted

Numerical evaluation of an elliptic integral in python

The problem is almost definitely with how QUADPACK (which is the backend used by scipy.integrate.quad) handles numerically small integrands. Essentially the ...
user avatar
  • 11.4k
8 votes
Accepted

Dynamically ending ODE integration in SciPy

The feature that you demand is called event location in Matlab ODE solvers pack, or rootfinding in SUNDIALS solvers suite terminology. Essentially this feature allows to stop integration exactly at ...
user avatar
8 votes
Accepted

Can numpy.linalg.solve use back substitution when possible?

Looking at the information of nympy.linalg.solve for dense matrices, it seems that they are calling LAPACK subroutine gesv, which perform the LU factorization of your matrix (without checking if the ...
user avatar
  • 409
8 votes

Plot integral function with scipy and matplotlib

First of all, your function $x\sin(\frac{1}{x})$ is singular in $x=0$. You might want to add an if clause like this: ...
user avatar
  • 6,046
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 ...
user avatar
  • 11.4k
8 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 ...
user avatar
8 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 ...
user avatar
  • 1,936
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 &...
user avatar
7 votes

Can numpy.linalg.solve use back substitution when possible?

No. The numpy.linalg.solve method uses LAPACK's DGESV, which is a general linear equation solver driver. If you know that your matrix is triangular, you should use a driver specialized for that matrix ...
user avatar
7 votes
Accepted

Solve non-linear set of three equations using scipy

Since your problem is small, you're probably best off trying fsolve or root. Both of these are interfaces to MINPACK and call <...
user avatar
7 votes
Accepted

Numerical solution of Geodesic differential equations with Python

The reason the resulting geodesic curve was deviating was because the calculated Christoffel symbol of second kind was incorrect. Using the correct Christoffel symbol : ...
user avatar
  • 425
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. ...
user avatar
  • 188
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 ...
user avatar
  • 383
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 ...
user avatar
  • 7,902
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 ...
user avatar
  • 3,789
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{...
user avatar
  • 876
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 ...
user avatar
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. <...
user avatar
  • 7,902
6 votes
Accepted

Solving a system of nonlinear PDEs by minimization

1) If you're just looking to solve the PDEs without any other optimization, then my answer would be "none of them". Algorithms that discretize partial differential equations and then solve them as ...
user avatar
6 votes

Scipy OdeInt solver with Neumann boundary conditions

From Ablowitz and Zeppetella we know that the analytical solution reads: \begin{equation} u(x,t)=\frac{1}{1+e^{{(-\frac{5}{6} t+ \frac{\sqrt{6}x}{6}})^2}} \end{equation} Usually, analytical ...
user avatar
6 votes

What does Python offer for distributed/parallel/GPU computing?

Here's a few options that are relatively easy to work with: One node - multiprocessing is the most straightforward thing to do. multiprocessing.map works well for ...
user avatar
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 ...
user avatar
  • 1,117
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 ...
user avatar
  • 11.4k
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 ...
user avatar

Only top scored, non community-wiki answers of a minimum length are eligible