Questions tagged [scipy]

SciPy is a Python-based ecosystem of open-source software for mathematics, science, and engineering.

Filter by
Sorted by
Tagged with
17
votes
1answer
5k views

When is Newton-Krylov not an appropriate solver?

Recently I have been comparing different non-linear solvers from scipy and was particularly impressed with the Newton-Krylov example in the Scipy Cookbook in which they solve a second order ...
14
votes
3answers
5k views

How are scientific computing workflows faring on Apple's M1 hardware

The initial wave of reviews for Apple's M1 hardware are out, and there's lots of generic benchmarks and data on workflows on professional programs for creative users, but I haven't seen anyone talking ...
13
votes
2answers
11k views

What is the fastest way to compute all eigenvalues of a very big and sparse adjacency matrix in python?

I'm trying to figure out if there is a faster way to compute all the eigenvalues and eigenvectors of a very big and sparse adjacency matrix than using scipy.sparse.linalg.eigsh As far as I know, this ...
13
votes
1answer
1k views

Can an approximated Jacobian with finite differences cause instability in the Newton method?

I have implemented a backward-Euler solver in python 3 (using numpy). For my own convenience and as an exercise, I also wrote a small function that computes a finite difference approximation of the ...
11
votes
5answers
10k views

Dynamically ending ODE integration in SciPy

I have a light ray moving through space-time, i.e. a curve in $\mathbb{R}^4$, parametrized by some variable λ. The exact trajectory, i.e. the coordinate functions $x^μ(λ)$ of the curve are given by ...
9
votes
1answer
627 views

Matrix Balancing Algorithm

I have been writing a control system toolbox from scratch and purely in Python3 (shameless plug : harold ). From my past research, I have always complaints about ...
8
votes
1answer
310 views

Forcing an ODE solver to preserve the norm

I have an ODE of the form $$ \frac{dy}{dt} = -i H y \enspace .$$ where $y$ is a complex vector and $H$ is a time dependent Hermitian matrix. The norm of the solution $y(t)$ at any point in time ...
8
votes
3answers
4k views

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

Using the SciPy/NumPy libraries, Python is a pretty cool and performing platform for scientific computing. I just wonder: When I have to go parallel (multi-thread, multi-core, multi-node, gpu), what ...
8
votes
1answer
4k views

Numerical evaluation of an elliptic integral in python

Goal: I need to evaluate numerically an integral of the following form: $$ \int_0^\infty \frac{dx}{(a^2+x)\sqrt{(a^2+x)(b^2+x)(c^2+x)}} $$ where $a,b,c \in \mathbb{R}$ are in the interval $(1,1000)$....
8
votes
1answer
2k views

Numerical solution of Geodesic differential equations with Python

I have made a solver based on the SymPy.diffgeom library, where I use Scipy.Integrate to solve the following system of second order differential equations : \begin{align} u'' &+ \Gamma^0_{00}(u')...
8
votes
2answers
1k views

Minimum path on known potential surface

I'm searching for the minimum path between the minima of a potential surface that is already known on a grid. (source: http://www.math.nus.edu.sg/~matrw/string/) Any point on the path is at an ...
7
votes
4answers
20k views

Plot integral function with scipy and matplotlib

I want to plot a numerical integral function of some function $f$ using scipy and matplotlib. How can I do this? I tried the ...
7
votes
4answers
6k views

Algorithm for high quality 1/f noise?

How can I generate arbitrarily high quality $1/f$ noise, for use in a model? My model involves a lot of feedback, over a large number of iterations, with a very high bandwidth, so I'd like the $1/f$ ...
7
votes
1answer
810 views

Why does LSODA fail to integrate the logistic function?

I'm comparing some of the different ODE integrators in scipy.integrate.ode on solving the logistic function: $$x(t) = \frac{1}{1+e^{-rt}}$$ $$\dot{x} = rx(1-x)$$ ...
6
votes
2answers
3k views

Is there an efficient way to form this block matrix with numpy or scipy?

Is there an efficient way to form this block matrix with numpy or scipy? $$ \left[ \begin{array}{cccc} \mathbf{B} & \mathbf{0} & \cdots & \mathbf{0}\\ \mathbf{AB} & \mathbf{B} & \...
6
votes
2answers
5k views

Solve non-linear set of three equations using scipy

I need to solve a non-linear set of three equations using scipy. However, I do not have any clue on which algorithm is suitable for my problem from a mathematical point of view (stability, ...
6
votes
2answers
2k views

Numerical Integration with Convergence Factor with SciPy: Problem with Improper Integral

I would like to perform the numerical integration of an integral of the form $$ \int_{-\infty}^\infty e^{i \omega 0+} G(i \omega, \mathbf{v}) d \omega ,$$ or, using the symmetry $G(i\omega)^* = G(-i \...
5
votes
3answers
297 views

Maximize a function of an orthogonal matrix

I'm trying to write up a small code that, given a set of normal vibrational modes for a molecule, will convert them to localized vibrational modes. To do this I'm following the procedure from J. Chem. ...
5
votes
1answer
163 views

Algorithm to factorize matrix whose many rows are already of upper triangular form?

I have a matrix whose many rows are already in the upper triangular form. $$\begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} & x_{5} \\ 0 & x_{22} & x_{23} & x_{...
5
votes
2answers
8k views

SciPy interpolation with Univariate Splines

I have coded a routine for interpolation with B-splines, only to discover later that this functionality is already included in Python's SciPy. However, I do not understand one parameter in the SciPy ...
5
votes
1answer
1k views

Solving for a set of coupled ODEs to get correct variable values

My question is about how I can solve a coupled system of ODE's, and print out the variables in a plot. I am solving for an q value and an e value, seen in this set of coupled ODE's below: $$ \begin{...
5
votes
5answers
9k views

scipy.sparse: Set row/column in sparse matrix to the identity without changing sparsity

I'm using the SciPy sparse.csr_matrix format for a finite element code. In applying the essential boundary conditions, I'm setting the desired value in the right ...
5
votes
1answer
1k views

Solving a system of nonlinear PDEs by minimization

I have two coupled nonlinear partial differential equations of the form: $ \begin{align} \dot{u} -f(u,u',u'',v,v',v'')=0 \\ \dot{v} -g(u,u',u'',v,v',v'')=0 \end{align} $ The boundary conditions are ...
5
votes
1answer
411 views

How to know which LAPACK's function is used by Scipy's eig function?

As far as I understood, scipy.linalg.eig use wrappers from scipy.lapack to compute the eigenvalues and eigenvectors of a matrix. ...
5
votes
1answer
97 views

How do I globally change the precision of a piece of code in Python to debug it?

I am solving a system of non-linear equations using the Newton-Raphson method in Python. This involves using the solve(Ax,b) function (...
5
votes
1answer
4k views

Solve rate equations with different reaction orders using SciPy ode

The image below illustrates the kinetic scheme I am trying to model. My first example focuses on the conversion of the wood-oil to non-volatiles and volatiles. The rate equations are as follows: $$ ...
5
votes
1answer
347 views

Quantify integration error of scipy ode / ODEPACK

I am trying to integrate a 2nd order ODE with potential several singularities using the lsoda solver wrapped in scipy.integrate.ode(). I would like to put an error bar on the solution or at least ...
5
votes
0answers
121 views

How to numerically evaluate this double Integral?

I want to evaluate the following integral: $$\int_{0}^{60} \ \left(\int_{0}^{2z} 0.5\cdot t \left(\mathrm{erf}(t-a) -1 \right)J_{0}(qt)\mathrm{d}t \right)^2 \mathrm{exp}\left(-\frac{(z-a)^2}{2s^2}\...
5
votes
0answers
400 views

Symmetric sparse direct solvers in scipy

scipy.linalg.solve, in its newer versions, has a parameter assume_a that can be used to specify that the matrix $A$ is symmetric ...
5
votes
0answers
426 views

Optimisation of matrix exponential

I have a 7000x7000 sparse matrix (scipy), which I want to exponentiate. I've tried using scipy.sparse.linalg.expm, which works quite well for smaller matrices (takes a few seconds for a 1000x1000 ...
4
votes
2answers
8k views

Tikhonov regularization in the non-negative least square - NNLS (python:scipy)

I am working on a project that I need to add a regularization into the NNLS algorithm. Is there a way to add the Tikhonov regularization into the NNLS implementation of scipy [1]? [2] talks about it, ...
4
votes
3answers
7k views

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

The question is if Python Numpy library can use back subsitution to solve Ax=b if possible, that is, if A is lower triangular? Do numerical linear algebra packages do this? I would think Numpy would ...
4
votes
2answers
333 views

Leveraging scipy for matrix free finite elements

This will be a very general question. I have a 3D finite element code in Python which I would like to extend to handle "large" problems (~10^8 unknowns in the global system). Right now I am using the ...
4
votes
4answers
10k views

LCM builtin in Python / Numpy

I can write a function to find LCM (lowest common multiple) of an array of integers, but I thought it must have been implemented in numpy or scipy and was expecting something like ...
4
votes
2answers
4k views

Number of equations and precision of SciPy's integrate.odeint()

Is there any reason why SciPy's integrate.odeint() should become less precise when the number of equations increase? I'm trying to solve these two sets of differential equations: $\frac{dy_1}{dx} = ...
4
votes
2answers
3k views

Numerical evaluation of a Gaussian Integral in Python?

Goal I'm trying to write code to compute the normalized Gaussian in the following, $$ \begin{equation} \int_{-\infty}^{\infty} \frac{1}{ \sigma \sqrt{2 \pi}} \exp\bigg( - \frac{(x - \mu)^{2}}{2 \...
4
votes
2answers
1k views

scipy odeint: sum of conservative ode equations does not remain zero as it is being solved — is this normal?

Let us say we have the following equations: dy1/dt = f(y1, t) [1] dy2/dt = g(y2, t) [2] The equations are such that they are "conservative", i.e. the ...
4
votes
2answers
385 views

Moore-Penrose pseudoinverse of singular rank degenerate matrix

I am trying to attain the Moore-Penrose pseudoinverse of a very large, very sparse, rank-degenerate, singular, and square matrix. ($75000 \times 75000$, near rank). The matrix is a graph Laplacian and ...
4
votes
2answers
16k views

scipy odeint - Excess work done on this call

I'm newbie both in calculus and Python/Scipy so I apologize if this question is too dumb. I'm trying to model flow between two pressure vessels. Let's say we have two points and a link between them ...
4
votes
1answer
96 views

Plotting optimum as a function of parameter in the objective

I am trying to minimize a 2d function using scipy.optimize. Specifically I want to plot the minimum value of the function fun as a function of the parameter wjk. The problem is that I cannot pass wjk ...
4
votes
1answer
193 views

Scipy Spline Interpolation Parameter

Documentation in scipy.interpolate (found at https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) states: "The parameter variable is given with the keyword argument, u, which ...
4
votes
0answers
394 views

Iteratively finding both left and right eigenvectors for non-symmetric complex matrix

I have a complex, non-Hermitian matrix $\mathbf{A}$, for which I need to find a few eigenvalues and eigenvectors in the generalised eigenvalue problem: $$\mathbf{A}\cdot \mathbf{x} = \lambda \mathbf{...
4
votes
0answers
79 views

Large residual when integrating 2nd order ode close to singularity with SciPy ode / ODEPACK

I am trying to integrate a 2nd order ODE with a singularity at close to the initial condition. Why do I get large residuals when I plug-in the result of my integration back into the ODE? The equation ...
3
votes
2answers
592 views

Computing numeric derivative via FFT - SciPy

I wrote the following code to compute the approximate derivative of a function using FFT: ...
3
votes
4answers
1k views

Numerical integration in Python with unknown constant

I’d like to solve the below equation for the unknown $T$: $$\int_0^\infty \frac{x^2}{\exp\left(\frac{x}{T}\right)-1}\kappa_x \mathrm{d}x = C,$$ where $C$ is a known constant and $\kappa_x$ is some ...
3
votes
2answers
2k views

Integration of the Fermi distribution using Python

I want to calculate the carrier concentration of my semiconductor using this equation: $$ n(x) = \frac{m^*}{\pi\hbar^2}\int_{E_k}^{\infty}\frac{1}{1+\exp\left(\frac{E-E_f}{k_BT}\right)} \mathrm{d}E $$...
3
votes
5answers
9k views

Fast c++ library to solve very big sparse systems

I am working on a project with electrical circuits, where I am trying to compute the voltages at all the nodes of an electrical circuit. I know that the electrical circuit is a perfect grid, so each ...
3
votes
2answers
2k views

Scipy OdeInt solver with Neumann boundary conditions

I'm using scipy.odeint to solve Fisher-Kolmogorov equation: \begin{equation} u_t = u_{xx}+u(1-u) \end{equation} The code can be found here. From Ablowitz and ...
3
votes
1answer
1k views

Numerical derivative in python

I am trying to take the numerical derivative of a dataset. My first attempt was to use the gradient function from numpy but in that case the graph of the derivative ...
3
votes
1answer
373 views

Get the roots of a Hermite interpolating polynomial

I am using Python 3.7 to write a program that requires me to calculate the root of the Hermite interpolating polynomial, given two points $\epsilon_0$, $\epsilon_1$, the function ($d(\epsilon_0)$,$d(\...