Speeding up python computation time (solving differential equations)

An alternative solution to the one of @IvanPerehiniak, is to use a JIT compiler like Numba so to do many low-level optimization that the CPython interpreter do not. Indeed, numerically intensive pure-Python code running on CPython are generally very inefficient. Numpy can provide relatively good performance for large arrays but it is very slow for small one. The thing is you use a lot of small arrays and pure-Python scalar operations. Numba is not a silver-bullet though: it just mitigate many overhead from Numpy and CPython. You still have to optimize the code further if you want to get a very fast code. Hopefully, this method can be combined with the one of @IvanPerehiniak (though the dictionary need not to be global which is cumbersome in many cases). Note Numba can pre-compute global constants for you. The compilation time is done during the first call or when the function has a user-defined explicit signature.

import numba as nb from scipy.integrate import quad from scipy.optimize import fsolve import numpy as np import matplotlib.pyplot as plt import scipy # These values can be changed masstot = 5.0 mass = 2.0 g= 9.8 l = 9.8 wan = (g/l)**(1/2) vuk = 0.1 oug = 1.0 @nb.njit def afad(lah): # Find first constant wan = 1.0 vuk = 0.1 oug = 1.0 kan = (12*(lah**4)*((3*(vuk**2)-(wan**2))))-((16*((wan**2)-(vuk**2))-(5*oug**2))*(lah**2))+(4*(oug**2)) return (kan) solua = fsolve(afad, 1) intsolua = np.sum(solua) @nb.njit def kfad(solua, wan, vuk): # Find second constant res = ((wan**2)-(vuk**2)-((2*(solua**2)*((2*(vuk**2))+(wan**2)))/((5*(solua**2))+4)))**(1/2) return (res) ksol = kfad(solua, wan, vuk) @nb.njit def deg(t, solua, vuk, ksol): # Find angle of pendulum relative to time res = 2*np.arctan(solua*np.exp(-1*vuk*t)*np.sin(ksol*t)) return(res) @nb.njit def chandeg(t, solua, vuk, ksol): # Find velocity of pendulum relative to time res = (((-2*solua*vuk*np.exp(vuk*t)*np.sin(ksol*t))+(2*solua*ksol*np.exp(vuk*t)*np.cos(ksol*t)))/(np.exp(2*vuk*t)+((solua**2)*(np.sin(ksol*t)**2)))) return(res) xs = np.linspace(0, 60, 20) # Value can be changed to alter plotting accuracy and length @nb.njit def dinte1(deg, bond, solua, vuk, ksol): # used to plot angle at at a certain time res = [] for x in (bond): res.append(deg(x, solua, vuk, ksol)) return res @nb.njit def dinte2(chandeg, bond, solua, vuk, ksol): # used to plot angular velocity at a certain time res = [] for x in (bond): res.append(chandeg(x, solua, vuk, ksol)) return res @nb.njit def dinte(a, bond, mass, l, solua, vuk, ksol, g, masstot ): # used to plot acceleration of system at certain time res = [] for x in (bond): res.append(a(x, mass, l, solua, vuk, ksol, g, masstot)) return res @nb.njit def a(t, mass, l, solua, vuk, ksol, g, masstot): # define acceleration of system to time return (((mass*l*(chandeg(t, solua, vuk, ksol)**2))+(mass*g*np.cos(deg(t, solua, vuk, ksol))))*np.sin(deg(t, solua, vuk, ksol))/masstot) # See: https://stackoverflow.com/questions/71244504/reducing-redundancy-for-calculating-large-number-of-integrals-numerically/71245570#71245570 @nb.cfunc('float64(float64)') def j(t): return np.sum(a(t, mass, l, intsolua, vuk, ksol, g, masstot)) j = scipy.LowLevelCallable(j.ctypes) # Cannot be jitted due to "quad" def f(ub): return quad(lambda ub: quad(j, 0, ub)[0], 0, ub)[0] # Cannot be jitted due to "f" not being jitted def int2(f, bond): # Integrates system acceleration twice to get posistion relative to time res = [] for x in (bond): res.append(f(x)) print(res) return res plt.plot(xs, int2(f, xs)) # This part of the program runs quite slowly #plt.plot(xs, dinte(a, xs, mass, l, solua, vuk, ksol, g, masstot)) #plt.plot(xs, dinte2(chandeg, xs, solua, vuk, ksol)) #plt.plot(xs, dinte1(deg, xs, solua, vuk, ksol)) plt.show() 

Here are results:

Initial solution: 35.5 s Ivan Perehiniak's solution: 5.9 s This solution (first run): 3.1 s This solution (second run): 1.5 s 

This solution is slower the first time the script is run because the JIT needs to compile all the functions the first time. Subsequent calls to the functions are significantly faster. In fact, int2 takes only 0.5 seconds on my machine the second time.