Questions tagged [scipy]

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

51 questions with no upvoted or accepted answers
Filter by
Sorted by
Tagged with
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
397 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
425 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
1answer
192 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
0answers
58 views

Does shift-invert method has invertibility issue?

Please note that I have nearly zero background on numerical methods. I understand the shift invert method as described in SciPy Tutorial The main argument of the above link is as follows. Suppose we ...
3
votes
0answers
55 views

Check if LinearOperator is symmetric

I have a scipy.sparse.linalg.LinearOperator object. I'd like to check if its associated matrix is symmetric without actually instantiating the matrix in the most ...
3
votes
0answers
48 views

Difference between wave vector and density matrix in numerical calculation of Schrödinger equation

I solved Schrödinger equation for a following tow-level time-dependent Hamiltonian numerically in two ways: ...
3
votes
0answers
76 views

Python routine to calculate shape resonances of H2

I am currently doing a project in which my aim is to write a program that can be used to calculate single and multi-channel shape resonances. So I'm looking at bound states and quasi-bound states. ...
3
votes
1answer
352 views

Poisson image blending artifacts

I am trying to implement Poisson image blending as in the paper Poisson Image Editing. This is the task of filling in a masked region of an image by minimizing $$\min_f\int_\Omega \left | \nabla f - \...
3
votes
0answers
260 views

Computing Small Eigenvalues with Sparse Symmetric Indefinite Mass Matrix

I want the eigenvalues of the following generalized eigenvalue problem: $$ Av = \lambda M v $$ where $A\in\mathbb{R}^{n\times n}$ is sparse, symmetric, and positive semi-definite $M\in\mathbb{R}^{n\...
3
votes
0answers
91 views

Inconsistency in optimize.minimize

I am trying to fit a time-dependent curve at each time step. I do so in minimizing along $x_c$ the quadratic error between the curve and a reference solution $ 1/(1 + \exp\left(\sqrt{S}(x-x_c)\right) $...
3
votes
0answers
197 views

Left eigenvectors using ARPACK

I'm trying to find both the dominant $k$ left and right eigenvectors, that is, $$V_L\mathcal{A} = \Lambda V_L\\ \mathcal{A}V_R = V_R\Lambda\\ V_LV_R = I_{k\times k}$$ $V_L$ being the $k\times N$ ...
3
votes
0answers
455 views

Gradients of non-uniformly sampled data in 3D space

I have measurements of magnetic field on a 3d grid. My measurements are distributed on four x-y planes similar to what is shown in the image below. The measurements roughly follow a Cartesian grid but ...
3
votes
0answers
2k views

BFGS Fails to Converge

The model I'm working on is a multinomial logit choice model. It's a very specific dataset so other existing MNLogit libraries don't fit with my data. So basically, it's a very complex function ...
2
votes
0answers
64 views

Providing Jacobian as LinearOperator in scipy.optimize.root

I asked this question a few days ago on stackoverflow, but I figure scicomp.stackexchange is probably a better place. Sorry for the double post. I want to solve a system of nonlinear equations using ...
2
votes
1answer
127 views

Good languages/packages for interior point optimization with non-linear constraints?

I'm currently using Python's scipy.optimize package to perform parameter estimation for a system of 10 ODEs. I have some observed data, and I'm trying to find the set of parameters which makes the ODE ...
2
votes
0answers
30 views

Scaling tensor approximation by symmetric tensor decomposition with SciPy's L-BFGS-B

I am trying to approximate a symmetric tensor of which the values are in the range of [1e-7,1e-4], by a symmetric tensor decomposition of lower rank. For this I am using the L-BFGS-B method in SciPy's ...
2
votes
0answers
57 views

Non-negative Least Squares to perform Inverse Laplace with weights

I'm trying to perform the inverse Laplace transform of a (noisy) dataset $y_i$ using Tikhonov regularization: $$\min \sum_{i=1}^{N} \left(\int_0^\infty e^{-s_i t} f(t) \, dt - y_i \right)^2 - \lambda^...
2
votes
0answers
42 views

ILUTP in sparse.linalg.spilu?

In Matlab, an ILU with threshold and pivoting (ILUTP) can be passed by default as: setup.type = 'ilutp'; [L, U] = ilu(A, setup); Looking for an equivalent in ...
2
votes
0answers
461 views

finding null space to a complex matrix

I need to solve the following equation: $$ \begin{pmatrix} \frac{\omega^2}{c^2}\varepsilon_x-\mu_z^{-1}k_y^2-\mu_y^{-1}k_z^2 & \mu_z^{-1}k_xk_y & \mu_y^{-1}k_xk_z\\ \mu_z^{-1}k_xk_y &\...
2
votes
0answers
285 views

How to reuse permutation-orderings within scipy's SuperLU-wrapper?

i'm solving sparse linear equations within scipy 0.18 which internally resorts to SuperLU (after umfpack got removed due to license-issues). Current, i'm doing a complete re-factorization in each ...
2
votes
0answers
176 views

Is it worth switching to timesteppers provided by PETSc if I can't write down a Jacobian for my problem? Case study with “the amoeba” toy problem

I am considering using petsc4py instead of scipy.integrate.odeint (which is a wrapper for Fortran solvers) for a problem ...
1
vote
1answer
170 views

Translating the Euler code in scipy's solve_ivp

My code is based on the similarity transformation X=VZ.I simulate the model for transformed equations involving Z by replacing ...
1
vote
0answers
51 views

Finite dimensional optimization problem over dynamical system

I am interested in solving numerically the following mathematical problem Consider an ode of the form $$ \dot q(t) = f(q(t),t_1,\ldots, t_N),\qquad t\in [0,T], $$ where $q\in \mathbb{R}^n$ is the ...
1
vote
0answers
1k views

Solve system of polynomial equations with Python

I have 5 at most 4th order polynomials in 5 variables, $$p_i(x_1,x_2,x_3,x_4,x_5) \qquad i = 1, \ldots, 5$$ where all coefficients are either rational or floating point. I'd would like to get the ...
1
vote
0answers
251 views

Rank filter on an nXm array using python

I would like to apply a rank filter on an nXm numpy array. Let's say I have this array: ...
1
vote
0answers
49 views

Update model parameter with new data, discarding old data

I have this dataset, and I am using y = (a * x^n) / (b + x^n) Hill function as the model, where a is the limit of the Hill curve,...
1
vote
0answers
339 views

Numerical integral of oscillating function with known zeros

I have a function that I need to numerically integrate from $0$ to $+\infty$, given by: $$I = \int_0^{+\infty} \mathrm{d}x\,x\,T^2(x)f(x)$$ where $T^2$ is an interpolated function that goes to $1$ ...
1
vote
0answers
107 views

solving numerically a 2D integral by using simps and quad combined

I have the following function $$f(x,y) =\left(\frac{1}{\exp(x-E_f)+1}-\frac{1}{\exp(x-y-E_f)+1} \right)\frac{1}{\sqrt{4-xy}}$$ I want to integrate it using Python 3. The domain of the $y$ variable ...
1
vote
0answers
667 views

Numerical Double integration with endpoint singularity in scipy Python gives incorrect answer

I am trying to integrate the following function in Python, $\int_{0}^{\infty}\int_{0}^{\infty} \dfrac{e^{-x-y}}{B(x,y)}dx dy$, where $B(x,y)$ is the beta function - $B(x,y) = \int_0^{1}a^{x-1}(1-a)^{...
1
vote
0answers
69 views

Convert scipy integration with one step to matlab integration

Scipy integration allows us to do ode integration one adaptive timestep at a time and do something to it. However, matlab ode needs us to specify a timespan , and determine the adaptive timestep ...
1
vote
0answers
151 views

vectorizating function integration scipy python

I need to implement the following in python: For a given discrete time series Zt (t=0 to T), find smallest t such that: $c\sum_{s=0}^t e^{[k(Zt-Zs)+m(t-s)]} >= \frac{p*}{1-p*} $ where c,k,m are ...
1
vote
0answers
901 views

linearly interpolate and determine gradients for data on non-uniform grid

I have measurements of a quantity on a 3d grid. My measurements are distributed on four x-y planes similar to what is shown in the image below. The measurements roughly follow a Cartesian grid but ...
1
vote
0answers
184 views

Generalized Hermite Function as eigenfunction of a differential operator

I'm going through this paper. The article defines function function $\phi_n^\mu(x)$ that is orthonormal on $L^2$ with measure $dm = dx$: \begin{equation} \phi^\mu_n =\left(\frac{\gamma_\mu(n)}{\...
1
vote
0answers
138 views

How to compute frank copula and its derivative accurately?

I need to fit a model using MLE with Frank copula by linking two discrete univate distribution function $u = F(x)$ and $v = F(y)$ together, and the joint distribution function is $$ \Phi(x,y) = C(F(x)...
0
votes
1answer
87 views

Error using scipy.integrate.solve_ivp: index error: the index 0 is out of bounds for axis 0 with size 0

I am recently working on the code based on the stick-slip phenomenon in Python. It's the stick-slip oscillator (chapter 3.4) in https://www.sciencedirect.com/science/article/pii/S0888327020301205. And ...
0
votes
0answers
60 views

L2 norm optimization problem

I have an optimization problem where i need to find an image x, that is very close to x' such that: monitor(x') is valid but monitor(x) is invalid. (output is valid when the neural network output is ...
0
votes
0answers
28 views

Set of integrators do not consistently solve an equation in Python

I must solve the following second order differential equation: $\delta \phi^{''}_{\mathbf{k}}+(3-\epsilon)\delta \phi^{'}_{\mathbf{k}}+\left(\frac{k^2}{a^2 H^2}+\frac{V_{,\phi\phi}}{H^2}-6\epsilon +4\...
0
votes
0answers
86 views

Applying the result of Cuthill-McKee in SciPy (followup)

This is a followup to Applying the result of Cuthill-McKee in SciPy , where I'm not sure the answer given is correct. It's also 4 years old, so I'm trying a new question. The question is still ...
0
votes
0answers
105 views

scipy's solve_ivp returns erroneous results for a stiff differential equation

I'm using scipy's solve_ivp for solving a stiff differential equation. I'm using method BDF for solving the same. I have already used MATLAB's ode23s and I'm getting correct results in MATLAB. However,...
0
votes
0answers
123 views

A parallelized GMRES solver?

My application calls for solving a dense, 40,000 x 40,000, ill-conditioned linear system. The native SciPy GMRES solver with preconditioning has worked well for my application and solving a single ...
0
votes
0answers
32 views

Comparing custom linear regression solver to SciPy equivalent in Python

From a given data set, I set out to complete a task which is below Fit the data of the previous exercise to fit Eq. (8.18) using the SciPy function ...
0
votes
0answers
46 views

Scipy.integrate.odeint is returning curves with almost the same frequency for different damping ratios, shouldn't they be different?

I am trying to solve the ODE for a harmonic oscillator using Scipy's odeint solver for different dampening factors. I'm using the following code, based off of this example: ...
0
votes
0answers
38 views

How to handle system of chemical reactions for a batch reactor SciPy solver

I have a system of chemical reactions where the rate equations represent a batch reactor model. The model is a system of ODEs which is solved with the SciPy ...
0
votes
0answers
43 views

Time sampling changes solution

I'm currently trying to solve a problem using numerical methods. The set-up is rather long, so I apologize in advance... TL;DR: My solutions change depending on how big my steps are and I don't know ...
0
votes
0answers
909 views

Tridiagonal Solver in Python

I have a code I'm working on that involves solving a 1D Schrodinger equation using a Crank-Nicolson time step. The code is written in NumPy/SciPy, and I was doing a bit of profiling and discovered ...
-1
votes
0answers
27 views

scipy.optimize.minimize(method=“SLSQP”) does not give the minimal value along its own search path

I would like to enquire about a situation I encountered while using the scipy function scipy.optimize.minimize(method="SLSQP"). It seems to loose track of the global minimum in the path of ...
-1
votes
0answers
24 views

Scipy: Sampling a custom distribution

Using Scipy, I would like to draw a random sample according to a certain probability density. The probability density function I consider is defined on the real line and depends on k parameters that I ...