I’ve spent a long time finally getting my algorithm in Python to invert Laplace transforms and live plot them to work, but unfortunately, although it seems to run fine within the function, it continuously just plots zeros!
I’ve wasted so much time getting this thing to work, and now I have no clue what could be wrong. So forgive me but I think I’ll have to post the whole thing.
The equation this is based on is a series-integral and it iterates until the series converges. However, it converges super slowly so I used Aitken’s delta-squared method to speed the whole thing up. The goal is to populate array y with each Aitken’s iteration from 0.000001 to 10 and then compare the estimation’ integral with the previous iteration’s if a random point in the former is close enough to the latter. This loops until the integrals have a difference of less than 0.01 or a set amount of iterations is reached (~30000).
This happens for 6 equations simultaneously using parallel processing. The graph live plots all 6 (in blue) compared to each of the actual equations they’re converging towards (in red) using matplotlib.
For debugging, I basically tried putting a bunch of print statements everywhere in the f_u function to see what was going on. Everything looked like it was going fine. But in the __main__ function and plot functions when I did the same, I was getting all zeros for the same array. I thought it was a local vs. global issue so I put global y in the plot and f_u, with no luck.
I know this isn’t doing this due to the parallelization; I tried without it with no luck. Basically, I suggest looking at what y does between the f_u function and the __main__ function, thats where the relationship with the y arrays breaks down I think, I just don’t know what to do to fix it.
To execute this, you’re gonna need all the packages in the import statement, including numba. You can remove the @jit decorations on the Laplace functions if you don’t want numba, but it will be sloooooow.
Here it is:
import matplotlib.pyplot as plt from matplotlib import style from matplotlib import animation import numpy as np from scipy.integrate import quad, simps from scipy.special import jv import math import random import datetime import multiprocessing as mp from numba import jit start = datetime.datetime.now() # Record time start # Initialize plot style.use('fivethirtyeight') fig, axes = plt.subplots(3, 2) ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = axes plt.autoscale(enable=True) for ax in axes.flat: ax.label_outer() # Initialize variables gamma = np.array([np.nan, 1, 1, 0, 1.5, 2, 1]) # Bromwich contour parameter for each equation num = 6 # Number of equations to be processed in total n = 0 # Repetition constant for equation 3, set to either zero or 1 (diverges at 1) 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 x_0 = 0 # Set first integral to 0 for now result = np.zeros(7) # Vector to store Aitken's iteration for closing() comparison previous = np.zeros(7) # Vector to store iteration before Aitken's for closing() comparison end = np.zeros(7) # Vector to store end conditions -- running = 0, ended = 1 u, step = np.linspace(0.000001, 10, max_count, retstep=True) # Initializing x values y = np.empty([num, 4, max_count]) # Initializing y values -- form is y[equation number][iteration in Aitken's][x value] y_denom = np.empty([num, max_count]) # Initializing array to store Aitken's denominator values for each equation x_iter = np.empty([num, max_count]) # Initializing array to store Aitken's iteration as starting point for next iteration def printer(num, k): # Prints the equation number, iteration and y value at 10 if num == 1: print(str(num) + " " + str(k) + " " + str(y[0][3][max_count-1]) + "\n") elif num == 2: print(str(num) + " " + str(k) + " " + str(y[1][3][max_count-1]) + "\n") elif num == 3: print(str(num) + " " + str(k) + " " + str(y[2][3][max_count-1]) + "\n") elif num == 4: print(str(num) + " " + str(k) + " " + str(y[3][3][max_count-1]) + "\n") elif num == 5: print(str(num) + " " + str(k) + " " + str(y[4][3][max_count-1]) + "\n") elif num == 6: print(str(num) + " " + str(k) + " " + str(y[5][3][max_count-1]) + "\n") def plot1(): # Plots estimation 1 (blue) and actual equation 1 (red) global y ax1.clear() line1, = ax1.plot(u, y[0][3][0:max_count]) ax1.plot(u, np.sin(u), 'tab:red') return line1, def plot2(): # Plots estimation 2 (blue) and actual equation 2 (red) global y ax2.clear() line2, = ax2.plot(u, y[1][3][0:max_count]) ax2.plot(u, np.heaviside(u, 1), 'tab:red') return line2, def plot3(): # Plots estimation 3 (blue) and actual equation 3 (red) global y ax3.clear() line3, = ax3.plot(u, y[2][3][0:max_count]) ax3.plot(u, np.cos(2 * np.sqrt(u)) / np.sqrt(np.pi * u), 'tab:red') return line3, def plot4(): # Plots estimation 4 (blue) and actual equation 4 (red) global y ax4.clear() line4, = ax4.plot(u, y[3][3][0:max_count]) ax4.plot(u, np.log(u), 'tab:red') return line4, def plot5(): # Plots estimation 5 (blue) and actual equation 5 (red) global y ax5.clear() line5, = ax5.plot(u, y[4][3][0:max_count]) ax5.plot(u, np.heaviside(u - 1, 1) - np.heaviside(u - 2, 1), 'tab:red') return line5, def plot6(): # Plots estimation 6 (blue) and actual equation 6 (red) global y ax6.clear() line6, = ax6.plot(u, y[5][3][0:max_count]) ax6.plot(u, jv(0, u), 'tab:red') return line6, def animate(i): # Plots all equations, used to live-update graph ln1, = plot1() ln2, = plot2() ln3, = plot3() ln4, = plot4() ln5, = plot5() ln6, = plot6() return ln1, ln2, ln3, ln4, ln5, ln6, @jit(nopython=True, cache=True) def f_p1(omega, u, k): # Equation 1, Laplace form b = (omega + k * np.pi) / u a = gamma[1] f1 = a * (1 / (a ** 2 + (b + 1) ** 2) + 1 / (a ** 2 + (b - 1) ** 2)) * np.cos(omega) return f1 @jit(nopython=True, cache=True) def f_p2(omega, u, k): # Equation 2, Laplace form b = (omega + k * np.pi) / u a = gamma[2] f2 = a / (a ** 2 + b ** 2) * np.cos(omega) return f2 @jit(nopython=True, cache=True) def f_p3(omega, u, k): # Equation 3, Laplace form b = (omega + k * np.pi) / u a = gamma[3] f3 = np.e ** (-a / (a ** 2 + b ** 2)) * (a ** 2 + b ** 2) ** (1 / 4) * np.cos( (math.atan2(b, a) + 2 * n * np.pi) / 2) * np.cos(b / (a ** 2 + b ** 2)) / ( np.sqrt(a ** 2 + b ** 2) * (np.cos((math.atan2(b, a) + 2 * n * np.pi) / 2) ** 2 + np.sin( (math.atan2(b, a) + 2 * n * np.pi) / 2))) * np.cos(omega) return f3 @jit(nopython=True, cache=True) def f_p4(omega, u, k): # Equation 4, Laplace form b = (omega + k * np.pi) / u a = gamma[4] f4 = -(a * np.log(a ** 2 + b ** 2) + a * np.e) / (a ** 2 + b ** 2) * np.cos(omega) return f4 @jit(nopython=True, cache=True) def f_p5(omega, u, k): # Equation 5, Laplace form b = (omega + k * np.pi) / u a = gamma[5] f5 = a * (np.e ** (-a) * np.cos(b) - np.e ** (-2 * a) * np.cos(2 * b)) / (a ** 2 + b ** 2) * np.cos(omega) return f5 @jit(nopython=True, cache=True) def f_p6(omega, u, k): # Equation 6, Laplace form b = (omega + k * np.pi) / u a = gamma[6] f6 = np.sqrt((a ** 2 - b ** 2 + 1) ** 2 + (2 * a * b) ** 2) * np.cos( (math.atan2(2 * a * b, a ** 2 - b ** 2 + 1) + 2 * n * np.pi) / 2) / ( ((a ** 2 - b ** 2 + 1) ** 2 + (2 * a * b) ** 2) * ( np.cos((math.atan2(2 * a * b, a ** 2 - b ** 2 + 1) + 2 * n * np.pi) / 2) ** 2 + np.sin( (math.atan2(2 * a * b, a ** 2 - b ** 2 + 1) + 2 * n * np.pi) / 2) ** 2)) * np.cos(omega) return f6 def closing(result, previous, num, 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 " + str(num) + " converges to " + str(diff) + " with " + str(k) + " iterations." + "\n") if num == 1: for i in range(max_count): y[0][3][i] += f_0(num, i) end[1] += 1 return 1 elif num == 2: for i in range(max_count): y[1][3][i] += f_0(num, i) end[2] += 1 return 1 elif num == 3: for i in range(max_count): y[2][3][i] += f_0(num, i) end[3] += 1 return 1 elif num == 4: for i in range(max_count): y[3][3][i] += f_0(num, i) end[4] += 1 return 1 elif num == 5: for i in range(max_count): y[4][3][i] += f_0(num, i) end[5] += 1 return 1 elif num == 6: for i in range(max_count): y[5][3][i] += f_0(num, i) end[6] += 1 return 1 else: print("Equation " + str(num) + " has not yet converged.\n") return 0 def f_0(num, i): # Initial integral in the series. Added at end after closing() convergence test if num == 1: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p1, 0, np.pi / 2, args=(u[i], 0))[0] elif num == 2: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p2, 0, np.pi / 2, args=(u[i], 0))[0] elif num == 3: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p3, 0, np.pi / 2, args=(u[i], 0))[0] elif num == 4: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p4, 0, np.pi / 2, args=(u[i], 0))[0] elif num == 5: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p5, 0, np.pi / 2, args=(u[i], 0))[0] elif num == 6: x_0 = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * quad( f_p6, 0, np.pi / 2, args=(u[i], 0))[0] return x_0 def f_u(u, num): # Main algorithm. Here, the summation integrals iterate global y, result, previous # Forces retrieval of from global of these arrays? (attempt to debug) k = 1 # Number of iterations for i in range(max_count): # First overall series iteration (initial 1st) if num == 1: y[0][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 2: y[1][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p2, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 3: y[2][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p3, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 4: y[3][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p4, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 5: y[4][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p5, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 6: y[5][0][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p6, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] for j in range(10000): # Iteration loop. Goes until closing() or k = ~30000 for each equation 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) if num == 1: y[0][0][i] = x_iter[0][i] elif num == 2: y[1][0][i] = x_iter[1][i] elif num == 3: y[2][0][i] = x_iter[2][i] elif num == 4: y[3][0][i] = x_iter[3][i] elif num == 5: y[4][0][i] = x_iter[4][i] elif num == 6: y[5][0][i] = x_iter[5][i] else: k += 1 for i in range(max_count): # Second iteration for Aitken's (2nd) if num == 1: y[0][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 2: y[1][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p2, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 3: y[2][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p3, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 4: y[3][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p4, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 5: y[4][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p5, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 6: y[5][1][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p6, -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) if num == 1: y[0][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p1, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 2: y[1][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p2, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 3: y[2][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p3, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 4: y[3][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p4, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 5: y[4][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p5, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] elif num == 6: y[5][2][i] = 2 * np.e ** (gamma[num] * u[i]) / (np.pi * u[i]) * (-1) ** k * quad( f_p6, -np.pi / 2, np.pi / 2, args=(u[i], k))[0] k += 1 for i in range(max_count): # Aitken's delta-squared method iteration using the previous 3 (4th) if num == 1: y_denom[0][i] = (y[0][2][i] - y[0][1][i]) - (y[0][1][i] - y[0][0][i]) if abs(y_denom[0][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[0][3][i] = y[0][2][i] - ((y[0][2][i] - y[0][1][i]) ** 2) / y_denom[0][i] elif num == 2: y_denom[1][i] = (y[1][2][i] - y[1][1][i]) - (y[1][1][i] - y[1][0][i]) if abs(y_denom[1][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[1][3][i] = y[1][2][i] - ((y[1][2][i] - y[1][1][i]) ** 2) / y_denom[1][i] elif num == 3: y_denom[2][i] = (y[2][2][i] - y[2][1][i]) - (y[2][1][i] - y[2][0][i]) if abs(y_denom[2][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[2][3][i] = y[2][2][i] - ((y[2][2][i] - y[2][1][i]) ** 2) / y_denom[2][i] elif num == 4: y_denom[3][i] = (y[3][2][i] - y[3][1][i]) - (y[3][1][i] - y[3][0][i]) if abs(y_denom[3][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[3][3][i] = y[3][2][i] - ((y[3][2][i] - y[3][1][i]) ** 2) / y_denom[3][i] elif num == 5: y_denom[4][i] = (y[4][2][i] - y[4][1][i]) - (y[4][1][i] - y[4][0][i]) if abs(y_denom[4][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[4][3][i] = y[4][2][i] - ((y[4][2][i] - y[4][1][i]) ** 2) / y_denom[4][i] elif num == 6: y_denom[5][i] = (y[5][2][i] - y[5][1][i]) - (y[5][1][i] - y[5][0][i]) if abs(y_denom[5][i]) < epsilon: print("Denominator too small to calculate. Exiting...\n") return 9 y[5][3][i] = y[5][2][i] - ((y[5][2][i] - y[5][1][i]) ** 2) / y_denom[5][i] k += 1 printer(num, k) 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 num == 1: if abs(y[0][3][rand] - y[0][2][rand]) < err: result[num] = simps(y[0][3][0:max_count], u, dx=step) #Integrates Aitken's iteration (4th) previous[num] = simps(y[0][2][0:max_count], u, dx=step) #Integrates integration before Aitken's (3rd) c1 = closing(result[num], previous[num], num, k) if c1 == 1: break elif num == 2: if abs(y[1][3][rand] - y[1][2][rand]) < err: result[num] = simps(y[1][3][0:max_count], u, dx=step) previous[num] = simps(y[1][2][0:max_count], u, dx=step) c2 = closing(result[num], previous[num], num, k) if c2 == 1: break elif num == 3: if abs(y[2][3][rand] - y[2][2][rand]) < err: result[num] = simps(y[2][3][0:max_count], u, dx=step) previous[num] = simps(y[2][2][0:max_count], u, dx=step) c3 = closing(result[num], previous[num], num, k) if c3 == 1: break elif num == 4: if abs(y[3][3][rand] - y[3][2][rand]) < err: result[num] = simps(y[3][3][0:max_count], u, dx=step) previous[num] = simps(y[3][2][0:max_count], u, dx=step) c4 = closing(result[num], previous[num], num, k) if c4 == 1: break elif num == 5: if abs(y[4][3][rand] - y[4][2][rand]) < err: result[num] = simps(y[4][3][0:max_count], u, dx=step) previous[num] = simps(y[4][2][0:max_count], u, dx=step) c5 = closing(result[num], previous[num], num, k) if c5 == 1: break elif num == 6: if abs(y[5][3][rand] - y[5][2][rand]) < err: result[num] = simps(y[5][3][0:max_count], u, dx=step) previous[num] = simps(y[5][2][0:max_count], u, dx=step) c6 = closing(result[num], previous[num], num, k) if c6 == 1: break for i in range(max_count): # Setting current Aitken's iteration to first iteration if num == 1: x_iter[0][i] = y[0][3][i] elif num == 2: x_iter[1][i] = y[1][3][i] elif num == 3: x_iter[2][i] = y[2][3][i] elif num == 4: x_iter[3][i] = y[3][3][i] elif num == 5: x_iter[4][i] = y[4][3][i] elif num == 6: x_iter[5][i] = y[5][3][i] # If iteration limit reached, sets process closing condition to 1 and returns to main end[num] = 1 return if __name__ == '__main__': i = 0 # Define and begin processes, one for each equation 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() # While any equation is running, update and show the plot every second while end[1] != 1 or end[2] != 1 or end[3] != 1 or end[4] != 1 or end[5] != 1 or end[6] != 1: animation.FuncAnimation(fig, animate, interval=1000, blit=True) plt.pause(1) plt.show(block=False) plt.pause(1) # Joining and closing each process when complete p1.join() p1.terminate() p2.join() p2.terminate() p3.join() p3.terminate() p4.join() p4.terminate() p5.join() p5.terminate() p6.join() p6.terminate() end = datetime.datetime.now() # Record time end print("Runtime: " + str(end - start)) # Updates and plots final result animate(i) plt.show() exit(0)
Hope you can help me out. Thanks in advance.
