1

我一直在研究一种用于在并行处理中反转拉普拉斯变换的算法（即同时集成 s 空间中的多个拉普拉斯函数），但我知道并行处理不是问题，因为这种方法为解决问题而实施。但是问题仍然存在，基本上我的代码一遍又一遍地吐出这个：

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)
4

1 回答 1

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)
于 2020-01-01T11:04:11.393 回答