在Python中并行化微分方程的分辨率

Question

我正在使用odeint 函数求解一个常微分方程组。是否有可能（如果是的话如何）轻松并行处理这类问题？

Answer 1

上面的答案是错误的，在数值上求解 ODE 需要每次迭代计算函数 f(t,y)=y' 多次，例如Runge-Kutta 四次。但我不知道任何 python 包这样做。

Answer 2

对 ODE 进行数值积分本质上是一种顺序运算，因为您需要每个结果来计算下一个结果（好吧，除非您从多个起点积分）。所以我猜答案是否定的。

Answer 3

编辑：哇，我刚刚意识到这个问题已经超过 3 年了。我仍然会留下我的答案，希望它能找到处于同样困境中的人的方式。很抱歉。

我遇到了同样的问题。我能够将这样的过程并行化如下。

首先你需要dispy。在那里，您会找到一些将为您执行并行化过程的程序。我不是dispy 方面的专家，但我使用它没有任何问题，而且我不需要配置任何东西。

那么，如何使用呢？

1. 运行python dispynode.py -d。如果您在运行主程序之前不运行此脚本，则不会执行并行作业。

2. 运行您的主程序。在这里，我发布了我使用的那个（对不起，一团糟）。您需要更改函数sim，并相应地更改您想要对结果执行的操作。不过，我希望我的程序可以作为您的参考。

   import os, sys, inspect
   
   #Add dispy to your path
   cmd_folder = os.path.realpath(os.path.abspath(os.path.split(inspect.getfile( inspect.currentframe() ))[0]))
   if cmd_folder not in sys.path:
       sys.path.insert(0, cmd_folder)
   
   # use this if you want to include modules from a subforder
   cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile( inspect.currentframe() ))[0],cmd_folder+"/dispy-3.10/")))
   if cmd_subfolder not in sys.path:
       sys.path.insert(0, cmd_subfolder)
   #----------------------------------------#
   #This function contains the differential equation to be simulated.    
   def sim(ic,e,O): #ic=initial conditions; e=Epsiolon; O=Omega 
       from scipy.integrate import ode
       import numpy as np
   
       #Diff Eq.
       def sys(t,x,e,O,z,b,l):
           p = 2.*e*O*np.sin(O*t)*(1-e*np.cos(O*t))/(z+(1-e*np.cos(O*t))**2)
           q = (1+4.*b/l*np.cos(O*t))*(z+(1-e*np.cos(O*t)))/( z+(1-e*np.cos(O*t))**2 )
           dx=np.zeros(2)
           dx[0] = x[1]
           dx[1] = -q*x[0]-p*x[1]
           return dx
       #Simulation.    
       t0=0; tEnd=10000.; dt=0.1
       r = ode(sys).set_integrator('dop853', nsteps=10,max_step=dt) #Definition of the integrator
       Y=[];S=[];T=[]
       # - parameters - # 
       z=0.5; l=1.0; b=0.06;
       # -------------- #
       color=1
       r.set_initial_value(ic, t0).set_f_params(e,O,z,b,l) #Set the parameters, the initial condition and the initial time
       #Loop to integrate.
       while r.successful() and r.t +dt < tEnd:
           r.integrate(r.t+dt)
           Y.append(r.y)
           T.append(r.t)
           if r.y[0]>1.25*ic[0]: #Bound. This is due to my own requirements.
               color=0
               break
           #r.y contains the solutions and r.t contains the time vector.
       return e,O,color #For each pair e,O return e,O and a color (0,1) which correspond to the color of the point in the stability chart (0=unstable) (1=stable)
       # ------------------------------------ #
   
   #MAIN PROGRAM where the parallel magic happens
   import matplotlib.pyplot as plt
   import dispy
   import numpy as np
   F=100 #Total files
   #Range of the values of Epsilon and Omega
   Epsilon = np.linspace(0,1,100)
   Omega_intervals   = np.linspace(0,4,F)
   
   ic=[0.1,0]
   
   cluster = dispy.JobCluster(sim) #This function sets that the cluster (array of processors) will be assigned the job sim.
   jobs = [] #Initialize the array of jobs
   
   for i in range(F-1):
       Data_Array=[]
       jobs = []
       Omega=np.linspace(Omega_intervals[i], Omega_intervals[i+1],10)
       print Omega
       for e in Epsilon:
           for O in Omega:
               job = cluster.submit(ic,e,O) #Send to the cluster a job with the specified parameters
               jobs.append(job) #Join all the jobs specified above
           cluster.wait()
       #Do the jobs
       for job in jobs:
           e,O,color = job()
           Data_Array.append([e,O,color])
   
       #Save the results of the simulation.
       file_name='Data'+str(i)+'.txt'
       f=open(file_name, 'a')
       f.write(str(Data_Array))
       f.close()