我是并行计算的新手,事实上,是数值方法的新手。我正在尝试使用以下形式的 python solve_ivp
求解微分方程:
y''(x) + (a^2 + x^2)y(x) = 0
y(0)=1
y'(0)=0
x=(0,100)
我想求解一系列 a
并将文件写入 a[i] y[i](80)
。
原始方程相当复杂,但本质上结构与上面定义的相同。我使用了 for
循环,它需要花费大量时间进行计算。在网上搜索我遇到了这个美丽的网站并找到了这个 question以及可能解决我面临的问题的相关答案。
我尝试了解决方案中提供的原始代码;但是,生成的输出未正确排序。 我的意思是,第二列的顺序不正确。
q1 a1 Y1
q1 a2 Y3
q1 a4 Y4
q1 a3 Y3
q1 a5 Y5
...
我什至尝试过使用一个参数的一个循环,但同样的问题仍然存在。下面是我的代码,使用相同的多处理方法,但使用 solve_ivp
import numpy as np
import scipy.integrate
import multiprocessing as mp
from scipy.integrate import solve_ivp
def fun(t, y):
# replace this function with whatever function you want to work with
# (this one is the example function from the scipy docs for odeint)
theta, omega = y
dydt = [omega, -a*omega - q*np.sin(theta)]
return dydt
#definitions of work thread and write thread functions
tspan = np.linspace(0, 10, 201)
def run_thread(input_queue, output_queue):
# run threads will pull tasks from the input_queue, push results into output_queue
while True:
try:
queueitem = input_queue.get(block = False)
if len(queueitem) == 3:
a, q, t = queueitem
sol = solve_ivp(fun, [tspan[0], tspan[-1]], [1, 0], method='RK45', t_eval=tspan)
F = 1 + sol.y[0].T[157]
output_queue.put((q, a, F))
except Exception as e:
print(str(e))
print("Queue exhausted, terminating")
break
def write_thread(queue):
# write thread will pull results from output_queue, write them to outputfile.txt
f1 = open("outputfile.txt", "w")
while True:
try:
queueitem = queue.get(block = False)
if queueitem[0] == "TERMINATE":
f1.close()
break
else:
q, a, F = queueitem
print("{} {} {} \n".format(q, a, F))
f1.write("{} {} {} \n".format(q, a, F))
except:
# necessary since it will throw an error whenever output_queue is empty
pass
# define time point sequence
t = np.linspace(0, 10, 201)
# prepare input and output Queues
mpM = mp.Manager()
input_queue = mpM.Queue()
output_queue = mpM.Queue()
# prepare tasks, collect them in input_queue
for q in np.linspace(0.0, 4.0, 100):
for a in np.linspace(-2.0, 7.0, 100):
# Your computations as commented here will now happen in run_threads as defined above and created below
# print('Solving for q = {}, a = {}'.format(q,a))
# sol1 = scipy.integrate.odeint(fun, [1, 0], t, args=( a, q))[..., 0]
# print(t[157])
# F = 1 + sol1[157]
input_tupel = (a, q, t)
input_queue.put(input_tupel)
# create threads
thread_number = mp.cpu_count()
procs_list = [mp.Process(target = run_thread , args = (input_queue, output_queue)) for i in range(thread_number)]
write_proc = mp.Process(target = write_thread, args = (output_queue,))
# start threads
for proc in procs_list:
proc.start()
write_proc.start()
# wait for run_threads to finish
for proc in procs_list:
proc.join()
# terminate write_thread
output_queue.put(("TERMINATE",))
write_proc.join()
请让我知道多处理中的错误,以便我可以在此过程中学习一些有关 python 中的多处理的知识。另外,如果有人让我知道在 python 中处理此类计算的最优雅/最有效的方法,我将不胜感激。谢谢
最佳答案
你想要的是一个 online 的东西排序。在这种情况下,您知道要显示结果的顺序(即,输入顺序),因此只需将输出累加到 priority queue 中即可。并在元素匹配您期望的下一个键时从中弹出元素。
一个没有显示并行性的小例子:
import heapq
def sort_follow(pairs,ref):
"""
Sort (a prefix of) pairs so as to have its first components
be the elements of (sorted) ref in order.
Uses memory proportional to the disorder in pairs.
"""
heap=[]
pairs=iter(pairs)
for r in ref:
while not heap or heap[0][0]!=r:
heapq.heappush(heap,next(pairs))
yield heapq.heappop(heap)
这种方法的优点——减少内存占用——对于只有几个 float 的小结果可能无关,但它很容易应用。
关于python - 在python中使用多处理编写微分方程输出的问题,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59728443/