如何使scipy.odeint的速度更快?

问题描述 投票:0回答:1

我目前正在解决一个559非线性微分方程的综合系统,我必须通过改变常量c1,c2 b和g来将解拟合到一些实验数据上。

代码是这样的。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import random as rd
from numba import jit

L=np.loadtxt('C:/Users/Pablo/Desktop/TFG/Probas/matriz_L_Pablo.txt')  
I=np.loadtxt('C:/Users/Pablo/Desktop/TFG/Probas/vector_I_Pablo.txt')



k=np.diag(L)
n=len(k)  #Contamos o numero de nodos

u=np.zeros(n)
for i in range (n):
    u[i]=rd.random()


M=np.zeros((n,n))
derivs=np.zeros(n)

c1=100 ; c2=10000 ; b=0.01 ; g=1


@jit
def f(y,t,params):
    suma=0
    c1,c2,b,g=params
    for i in range(n):
        for j in range(n):
            if i==j:
                M[i,i]=(1-y[i]/b)+g*(1-y[i])+c2*I[i]*(1/n-1)
            if i!=j:
                M[i,j]=(1/n)*(c1*L[i,j]+c2*I[i])
            out=(M[i,j]*y[j])
            suma=suma+out
        derivs[i]=suma
        suma=0
    return derivs





#Condicions iniciais

y0=u

#lista cos parametros
params=[c1,c2,b,g]

#tempos de int
tf=1
deltat=0.001
t=np.arange(0,tf,deltat)

#solucion

sol= odeint(f, y0,t, args=(params,))

(对不起,如果不是很清楚,我是第一次来这里)

optimization scipy odeint
1个回答
0
投票

你可以尝试将你的代码矢量化。函数f做了2件事情--首先它创建矩阵M,然后做乘法$$M。y$$. 乘法$$My$$很容易向量化,因为我们只需要使用numpy的matmul函数。

def f(y,t,params):
    suma=0
    c1,c2,b,g=params
    for i in range(n):
        for j in range(n):
            if i==j:
                M[i,i]=(1-y[i]/b)+g*(1-y[i])+c2*I[i]*(1/n-1)
            if i!=j:
                M[i,j]=(1/n)*(c1*L[i,j]+c2*I[i])

    return np.matmul(M,y)

这对运行时间应该有一定的帮助。但最耗时的部分是,每次调用f时,整个矩阵M都会形成,而且每次只形成一个元素。在调用f时,M中唯一需要修改的部分,是依赖于y的部分,所以M中所有的偏对角线项都可以在调用ODE求解器之前被填充。因此,如果M是(569x569),你不必每次调用f时都要计算M的所有~250000+元素,而只需要计算M中对角线上的569个元素。M 的其余 250000 个元素不依赖于 y,并且在调用 Ode 解算器之前指定。做这个修改应该会带来巨大的速度提升,因为这似乎是你代码中的主要瓶颈。最后,你也可以通过使用类似 numpy.diag_indices 这样的方法对 M 的对角线进行向量填充。

© www.soinside.com 2019 - 2024. All rights reserved.