将 linspace 向量发送到函数会使向量在函数开始之前全为零
Sending a linspace vector to a function makes the vector all zeroes before the function starts
我一直在研究一种在并行处理中反转拉普拉斯变换的算法（即同时集成 s-space 中的多个拉普拉斯函数），但我知道事实上，并行处理不是问题，因为此方法是为了解决问题而实施的 。但问题仍然存在，基本上我的代码一遍又一遍地吐出这个：
Denominator too small. Exiting...
以及各种域错误。我指出了问题所在。基本上是这样的：
for u in np.linspace(0.000001, 100, 1000000):
在 u 传递给我想要集成的任何函数时变为全零。不确定是不是积分功能：
x_1[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
或拉普拉斯函数 f_p 本身，但我已经尝试了一切，我不能浪费更多时间来旋转我的轮子。如果能得到任何帮助，我将不胜感激。这是 运行 第一个拉普拉斯反演所需的代码：
import numpy as np
from scipy.integrate import quad
import math
import multiprocessing as mp
gamma = 10
num = 1
k = 1
n = 0
epsilon = 10 ** -16
max_err = 10 ** -10
max_count = 1000000
x_iter = np.arange(1, 7)
x_1 = np.arange(1, 7)
x_2 = np.arange(1, 7)
x_3 = np.arange(1, 7)
denom = np.arange(1, 7)
AX = np.arange(1, 7)
a = gamma
def f_p1(omega, u):
b = (omega + k * math.pi) / u
f1 = a * (1 / (a ** 2 + (b + 1) ** 2) + 1 / (a ** 2 + (b - 1) ** 2))
return f1 * math.cos(omega)
def f_p(num):
if num == 1:
return f_p1
def f_0(u, num):
x = 2 * math.e ** (gamma * u) / (math.pi * u)
return x * quad(f_p(num), 0, math.pi / 2, args=(u,))[0]
def f_u(u, num):
k = 1
x = 2 * math.e ** (gamma * u) / (math.pi * u)
x_1[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
AX[num] = f_0(u, num)
for j in range(max_count):
if j > 0:
x_1[num] = x_iter[num]
k += 1
x_2[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
k += 1
x_3[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
denom[num] = (x_3[num] - x_2[num]) - (x_2[num] - x_1[num])
if(abs(denom[num]) < epsilon):
print("Denominator too small. Exiting...\n")
break
AX[num] = x_3[num] - ((x_3[num] - x_2[num]) ** 2) / denom[num]
k += 1
if(abs(AX[num] - x_3[num]) < max_err):
print("Equation " + str(num) + "converges to " + str(AX[num]) + "\n")
break
x_iter[num] = AX[num]
np.zeros(x_1)
np.zeros(x_2)
np.zeros(x_3)
np.zeros(denom)
np.zeros(AX)
if __name__ == '__main__':
mp.set_start_method('spawn')
for u in np.linspace(0.000001, 100, max_count):
p1 = mp.Process(target=f_u, args=(u, 1))
p2 = mp.Process(target=f_u, args=(u, 2))
p3 = mp.Process(target=f_u, args=(u, 3))
p4 = mp.Process(target=f_u, args=(u, 4))
p5 = mp.Process(target=f_u, args=(u, 5))
p6 = mp.Process(target=f_u, args=(u, 6))
p1.start()
p2.start()
p3.start()
p4.start()
p5.start()
p6.start()
p1.join()
p2.join()
p3.join()
p4.join()
p5.join()
p6.join()
exit(0)
这是一个愚蠢的问题，我没有真正理解 scipy 的 quad 函数是如何工作的，不小心把 f_1(num) 放在了我想要函数集成的地方关于 omega。我也没有将 k 放入参数中，而是在我的 Aitken 加速后添加而不是从 k 中减去（等于 AX 的部分）。
最后我也有np.zeros(AX)，这也导致了一堆问题。总的来说就是乱七八糟。
它现在可以工作了，我按照那个人的建议删除了多处理，以便同时绘制这些函数。没有看到性能上有太大的影响，去看看吧。
但是，出于某种原因它仍然没有给出正确答案，我知道这不是拉普拉斯方程（f_1）不正确。我认为这与我遍历向量 u 的方式有关... 如果有人知道如何针对变量进行积分并在 args= 部分中有另一个变量那个没有集成，让我知道！
这是它现在的样例，供后代使用。最初的方程式是为了对多个方程式执行此操作而构建的，但我将其减少为 1 以使其相对简短。也尝试使用 gamma 参数——它改变了收敛的方式，但不是我想要的方式：
import matplotlib.pyplot as plt
from matplotlib import style
import numpy as np
from scipy.integrate import quad, simps
import random
import datetime
from numba import jit
start = datetime.datetime.now() # Record time start
# Initialize variables
gamma = 2 # Bromwich contour parameter
max_count = 1000 # Number of data points in x and y
epsilon = 10 ** -40 # Minimum number to divide by in Aitken's iteration
err = 10 ** -4 # Maximum difference in y to trigger comparison in closing()
max_err = 0.01 # Maximum difference in integrals computed in closing() to end iteration process
result = 0 # Store Aitken's iteration for closing() comparison
previous = 0 # Store iteration before Aitken's for closing() comparison
u, step = np.linspace(0.000005, 10, max_count, retstep=True) # Initializing x values
y = np.empty([4, max_count]) # Initializing y values -- form is y[iteration in Aitken's][x value]
y_denom = np.empty([max_count]) # Initializing array to store Aitken's denominator values
x_iter = np.empty([max_count]) # Initializing array to store Aitken's iteration as starting point for next iteration
# Initialize plot
style.use('fivethirtyeight')
fig = plt.figure()
ax = fig.add_subplot(111)
plt.autoscale(enable=True)
def printer(k): # Prints the iteration and y value at x = 10
print(str(k) + " " + str(y[3][max_count-1]) + "\n")
def plot1(): # Plots estimation 1 (blue) and actual equation 1 (red)
ax.clear()
line1, = ax.plot(u, y[3][:])
ax.plot(u, np.sin(u), 'tab:red')
return line1,
@jit(nopython=True, cache=True)
def f_p1(omega, u, k): # Equation 1, Laplace form
b = (omega + k * np.pi) / u
a = gamma
f1 = a * (1 / (a ** 2 + (b + 1) ** 2) + 1 / (a ** 2 + (b - 1) ** 2)) * np.cos(omega)
return f1
def closing(result, previous, k): # Checks for convergence. If so, adds each term to f_0 and ends process
if abs(result - previous) < max_err:
diff = abs(result - previous)
print("Equation converges to " + str(diff) + " with " + str(k) + " iterations." + "\n")
for i in range(max_count):
y[3][i] += f_0(i)
return 1
else:
print("Equation has not yet converged.\n")
return 0
def f_0(i): # Initial integral in the series. Added at end after closing() convergence test
x_0 = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * quad(
f_p1, 0, np.pi / 2, args=(u[i], 0))[0]
return x_0
def f_u(u): # Main algorithm. Here, the summation integrals iterate
k = 1 # Number of iterations
for i in range(max_count): # First overall series iteration (initial 1st)
y[0][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
for j in range(2000): # Iteration loop. Goes until closing() or k = ~2000 (~8000 total iterations)
if j > 0:
for i in range(max_count): # Sets the first iteration (out of 3 for Aitken's) equal to the previous Aitken's or the first iteration if j = 0 (1st)
y[0][i] = x_iter[i]
else:
k += 1
for i in range(max_count): # Second iteration for Aitken's (2nd)
y[1][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
k += 1
for i in range(max_count): # Third iteration for Aitken's (3rd)
y[2][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
for i in range(max_count): # Aitken's delta-squared method iteration using the previous 3 (4th)
y_denom[i] = (y[2][i] - y[1][i]) - (y[1][i] - y[0][i])
if abs(y_denom[i]) < epsilon:
print("Denominator too small to calculate. Exiting...\n")
return 9
y[3][i] = y[2][i] - ((y[2][i] - y[1][i]) ** 2) / y_denom[i]
k -= 1
printer(k)
plot1()
rand = random.randrange(0, max_count) # Setting random number for x value
# Comparison between y at random x in Aitken's vs the previous iteration. If close enough, triggers closing()
if abs(y[3][rand] - y[2][rand]) < err:
result = simps(y[3][:], u, dx=step) #Integrates Aitken's iteration (4th)
previous = simps(y[2][:], u, dx=step) #Integrates integration before Aitken's (3rd)
c1 = closing(result, previous, k)
if c1 == 1:
return 0
for i in range(max_count): # Setting current Aitken's iteration to first iteration
x_iter[i] = y[3][i]
return
if __name__ == '__main__':
# Run the algorithm and plot
f_u(u)
plot1()
end = datetime.datetime.now() # Record time end
print("Runtime: " + str(end - start))
# Plots final result
plt.show()
exit(0)
我一直在研究一种在并行处理中反转拉普拉斯变换的算法（即同时集成 s-space 中的多个拉普拉斯函数），但我知道事实上，并行处理不是问题，因为此方法是为了解决问题而实施的 。但问题仍然存在，基本上我的代码一遍又一遍地吐出这个：
Denominator too small. Exiting...
以及各种域错误。我指出了问题所在。基本上是这样的：
for u in np.linspace(0.000001, 100, 1000000):
在 u 传递给我想要集成的任何函数时变为全零。不确定是不是积分功能：
x_1[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
或拉普拉斯函数 f_p 本身，但我已经尝试了一切，我不能浪费更多时间来旋转我的轮子。如果能得到任何帮助，我将不胜感激。这是 运行 第一个拉普拉斯反演所需的代码：
import numpy as np
from scipy.integrate import quad
import math
import multiprocessing as mp
gamma = 10
num = 1
k = 1
n = 0
epsilon = 10 ** -16
max_err = 10 ** -10
max_count = 1000000
x_iter = np.arange(1, 7)
x_1 = np.arange(1, 7)
x_2 = np.arange(1, 7)
x_3 = np.arange(1, 7)
denom = np.arange(1, 7)
AX = np.arange(1, 7)
a = gamma
def f_p1(omega, u):
b = (omega + k * math.pi) / u
f1 = a * (1 / (a ** 2 + (b + 1) ** 2) + 1 / (a ** 2 + (b - 1) ** 2))
return f1 * math.cos(omega)
def f_p(num):
if num == 1:
return f_p1
def f_0(u, num):
x = 2 * math.e ** (gamma * u) / (math.pi * u)
return x * quad(f_p(num), 0, math.pi / 2, args=(u,))[0]
def f_u(u, num):
k = 1
x = 2 * math.e ** (gamma * u) / (math.pi * u)
x_1[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
AX[num] = f_0(u, num)
for j in range(max_count):
if j > 0:
x_1[num] = x_iter[num]
k += 1
x_2[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
k += 1
x_3[num] = x * (-1) ** k * quad(f_p(num), -math.pi / 2, math.pi / 2, args=(u,))[0]
denom[num] = (x_3[num] - x_2[num]) - (x_2[num] - x_1[num])
if(abs(denom[num]) < epsilon):
print("Denominator too small. Exiting...\n")
break
AX[num] = x_3[num] - ((x_3[num] - x_2[num]) ** 2) / denom[num]
k += 1
if(abs(AX[num] - x_3[num]) < max_err):
print("Equation " + str(num) + "converges to " + str(AX[num]) + "\n")
break
x_iter[num] = AX[num]
np.zeros(x_1)
np.zeros(x_2)
np.zeros(x_3)
np.zeros(denom)
np.zeros(AX)
if __name__ == '__main__':
mp.set_start_method('spawn')
for u in np.linspace(0.000001, 100, max_count):
p1 = mp.Process(target=f_u, args=(u, 1))
p2 = mp.Process(target=f_u, args=(u, 2))
p3 = mp.Process(target=f_u, args=(u, 3))
p4 = mp.Process(target=f_u, args=(u, 4))
p5 = mp.Process(target=f_u, args=(u, 5))
p6 = mp.Process(target=f_u, args=(u, 6))
p1.start()
p2.start()
p3.start()
p4.start()
p5.start()
p6.start()
p1.join()
p2.join()
p3.join()
p4.join()
p5.join()
p6.join()
exit(0)
这是一个愚蠢的问题，我没有真正理解 scipy 的 quad 函数是如何工作的，不小心把 f_1(num) 放在了我想要函数集成的地方关于 omega。我也没有将 k 放入参数中，而是在我的 Aitken 加速后添加而不是从 k 中减去（等于 AX 的部分）。
最后我也有np.zeros(AX)，这也导致了一堆问题。总的来说就是乱七八糟。
它现在可以工作了，我按照那个人的建议删除了多处理，以便同时绘制这些函数。没有看到性能上有太大的影响，去看看吧。
但是，出于某种原因它仍然没有给出正确答案，我知道这不是拉普拉斯方程（f_1）不正确。我认为这与我遍历向量 u 的方式有关... 如果有人知道如何针对变量进行积分并在 args= 部分中有另一个变量那个没有集成，让我知道！
这是它现在的样例，供后代使用。最初的方程式是为了对多个方程式执行此操作而构建的，但我将其减少为 1 以使其相对简短。也尝试使用 gamma 参数——它改变了收敛的方式，但不是我想要的方式：
import matplotlib.pyplot as plt
from matplotlib import style
import numpy as np
from scipy.integrate import quad, simps
import random
import datetime
from numba import jit
start = datetime.datetime.now() # Record time start
# Initialize variables
gamma = 2 # Bromwich contour parameter
max_count = 1000 # Number of data points in x and y
epsilon = 10 ** -40 # Minimum number to divide by in Aitken's iteration
err = 10 ** -4 # Maximum difference in y to trigger comparison in closing()
max_err = 0.01 # Maximum difference in integrals computed in closing() to end iteration process
result = 0 # Store Aitken's iteration for closing() comparison
previous = 0 # Store iteration before Aitken's for closing() comparison
u, step = np.linspace(0.000005, 10, max_count, retstep=True) # Initializing x values
y = np.empty([4, max_count]) # Initializing y values -- form is y[iteration in Aitken's][x value]
y_denom = np.empty([max_count]) # Initializing array to store Aitken's denominator values
x_iter = np.empty([max_count]) # Initializing array to store Aitken's iteration as starting point for next iteration
# Initialize plot
style.use('fivethirtyeight')
fig = plt.figure()
ax = fig.add_subplot(111)
plt.autoscale(enable=True)
def printer(k): # Prints the iteration and y value at x = 10
print(str(k) + " " + str(y[3][max_count-1]) + "\n")
def plot1(): # Plots estimation 1 (blue) and actual equation 1 (red)
ax.clear()
line1, = ax.plot(u, y[3][:])
ax.plot(u, np.sin(u), 'tab:red')
return line1,
@jit(nopython=True, cache=True)
def f_p1(omega, u, k): # Equation 1, Laplace form
b = (omega + k * np.pi) / u
a = gamma
f1 = a * (1 / (a ** 2 + (b + 1) ** 2) + 1 / (a ** 2 + (b - 1) ** 2)) * np.cos(omega)
return f1
def closing(result, previous, k): # Checks for convergence. If so, adds each term to f_0 and ends process
if abs(result - previous) < max_err:
diff = abs(result - previous)
print("Equation converges to " + str(diff) + " with " + str(k) + " iterations." + "\n")
for i in range(max_count):
y[3][i] += f_0(i)
return 1
else:
print("Equation has not yet converged.\n")
return 0
def f_0(i): # Initial integral in the series. Added at end after closing() convergence test
x_0 = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * quad(
f_p1, 0, np.pi / 2, args=(u[i], 0))[0]
return x_0
def f_u(u): # Main algorithm. Here, the summation integrals iterate
k = 1 # Number of iterations
for i in range(max_count): # First overall series iteration (initial 1st)
y[0][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
for j in range(2000): # Iteration loop. Goes until closing() or k = ~2000 (~8000 total iterations)
if j > 0:
for i in range(max_count): # Sets the first iteration (out of 3 for Aitken's) equal to the previous Aitken's or the first iteration if j = 0 (1st)
y[0][i] = x_iter[i]
else:
k += 1
for i in range(max_count): # Second iteration for Aitken's (2nd)
y[1][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
k += 1
for i in range(max_count): # Third iteration for Aitken's (3rd)
y[2][i] = 2 * np.e ** (gamma * u[i]) / (np.pi * u[i]) * (-1) ** k * quad(
f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0]
for i in range(max_count): # Aitken's delta-squared method iteration using the previous 3 (4th)
y_denom[i] = (y[2][i] - y[1][i]) - (y[1][i] - y[0][i])
if abs(y_denom[i]) < epsilon:
print("Denominator too small to calculate. Exiting...\n")
return 9
y[3][i] = y[2][i] - ((y[2][i] - y[1][i]) ** 2) / y_denom[i]
k -= 1
printer(k)
plot1()
rand = random.randrange(0, max_count) # Setting random number for x value
# Comparison between y at random x in Aitken's vs the previous iteration. If close enough, triggers closing()
if abs(y[3][rand] - y[2][rand]) < err:
result = simps(y[3][:], u, dx=step) #Integrates Aitken's iteration (4th)
previous = simps(y[2][:], u, dx=step) #Integrates integration before Aitken's (3rd)
c1 = closing(result, previous, k)
if c1 == 1:
return 0
for i in range(max_count): # Setting current Aitken's iteration to first iteration
x_iter[i] = y[3][i]
return
if __name__ == '__main__':
# Run the algorithm and plot
f_u(u)
plot1()
end = datetime.datetime.now() # Record time end
print("Runtime: " + str(end - start))
# Plots final result
plt.show()
exit(0)