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)...
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. ...
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:
...
Only top voted, non community-wiki answers of a minimum length are eligible
