New answers tagged

5

You could define a linear opearator and pass it to the function eigsh. Ideally, your matrices $L$ and $X$ are sparse so you can take advantage of the matrix-vector product. In your case, you would have something like the following. import numpy as np from scipy.sparse.linalg import LinearOperator, eigsh def mv(v): a = 2.3 return L@(L@v) + a*X@(X.T@v)...


0

Rejection sampling is straightforward to implement for this case. import numpy as np def rng_xpy(n, rng=None, chunk_size=1024): rng = np.random.default_rng(rng) rvs = [] n_drawn = 0 while n_drawn < n: # Draw numbers from U(0, 1) for x, y, and z # We draw numbers in chunks for efficiency x, y, z = rng.random((3, ...


2

All-over this is a nicely structured code. The main problems are related to the Runge-Kutta solvers, where the first-order system was not uniformly applied to the computation. What is obviously wrong can be found in the first lines of the solver file def rk2_derivatives(edo, qk, pk, dt, bodies): k1 = dt * edo(qk, pk, bodies) k2 = dt * edo(qk + (dt * k1), ...


1

import numpy as np def vectors(max_n,max_t,k=0.05,l=0.01): times = np.arange(1,max_t,1) ##exclude t=0, set in initial condition indicies = np.arange(0,max_n,1) #initial conditions: F = np.zeros((max_n,max_t)) F[0,0] = 1 for t in times: for n in indicies: F[n,t] = -(l+(k*n)-1)*F[n,t-1] if n!=max_n-1: ...


3

Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. ...


Top 50 recent answers are included