I am using the trapezoidal method to integrate the squared modulus of this Gaussian function:
In fact, the squared modulus of this function can be integrated analytically, resulting in:
When I try to solve the integral both analytically and numerically, I get the following result (I simplify by setting ( t_0 = 0 ) and ( \sigma_t = 1 )):
Essentially, there is an offset on the y-axis. This concerns me because I then have to subtract the obtained values from another function, and this offset propagates throughout the code.
I can safely use the analytical solution for one case, but in the other, since the function comes from a frequency-domain reflection relationship of this Gaussian pulse, I am afraid that even in that case, when I integrate the squared modulus, there might be an offset as well.
Do you know what could be causing this and could you possibly suggest what to do? This is my first time dealing with numerical methods.
from scipy.integrate import trapezoid
from scipy.interpolate import CubicSpline
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
def u_time(t_0):
prefactor = 1 / (np.pi**0.25)
def func(time):
return prefactor * np.exp(-(time - t_0) ** 2 / 2)
return func
def u_integral(t_0, times):
u_values = u_time(t_0)(times)
u_abs_2 = np.abs(u_values) ** 2
u_abs_2_int = np.zeros_like(times, dtype=npplex128)
for k in range(1, len(times)):
u_abs_2_int[k] = trapezoid(u_abs_2[:k], times[:k])
return CubicSpline(times, u_abs_2_int)
def gaussian_squared_integral(t_0):
return lambda time: 0.5 * (erf(time - t_0) + erf(t_0))
t_0 = 0
time_arr = np.linspace(-20, 100, 1000)
analytic = u_integral(t_0, time_arr)(time_arr)
numeric = gaussian_squared_integral(t_0)(time_arr)
plt.plot(time_arr, analytic, label="Analytic")
plt.plot(time_arr, numeric, label="Numeric")
plt.legend()
plt.show()
I am using the trapezoidal method to integrate the squared modulus of this Gaussian function:
In fact, the squared modulus of this function can be integrated analytically, resulting in:
When I try to solve the integral both analytically and numerically, I get the following result (I simplify by setting ( t_0 = 0 ) and ( \sigma_t = 1 )):
Essentially, there is an offset on the y-axis. This concerns me because I then have to subtract the obtained values from another function, and this offset propagates throughout the code.
I can safely use the analytical solution for one case, but in the other, since the function comes from a frequency-domain reflection relationship of this Gaussian pulse, I am afraid that even in that case, when I integrate the squared modulus, there might be an offset as well.
Do you know what could be causing this and could you possibly suggest what to do? This is my first time dealing with numerical methods.
from scipy.integrate import trapezoid
from scipy.interpolate import CubicSpline
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
def u_time(t_0):
prefactor = 1 / (np.pi**0.25)
def func(time):
return prefactor * np.exp(-(time - t_0) ** 2 / 2)
return func
def u_integral(t_0, times):
u_values = u_time(t_0)(times)
u_abs_2 = np.abs(u_values) ** 2
u_abs_2_int = np.zeros_like(times, dtype=npplex128)
for k in range(1, len(times)):
u_abs_2_int[k] = trapezoid(u_abs_2[:k], times[:k])
return CubicSpline(times, u_abs_2_int)
def gaussian_squared_integral(t_0):
return lambda time: 0.5 * (erf(time - t_0) + erf(t_0))
t_0 = 0
time_arr = np.linspace(-20, 100, 1000)
analytic = u_integral(t_0, time_arr)(time_arr)
numeric = gaussian_squared_integral(t_0)(time_arr)
plt.plot(time_arr, analytic, label="Analytic")
plt.plot(time_arr, numeric, label="Numeric")
plt.legend()
plt.show()
- What time are you integrating from? One of them is the integral from -50; the other is the integral from 0. – lastchance Commented Mar 21 at 15:29
- It was my bad, I just posted the wrong code! – Frostman Commented Mar 21 at 16:12
1 Answer
Reset to default 3In fact, the squared modulus of this function can be integrated analytically, resulting in:
It's that plus C. When you take the indefinite integral of this, you get that value, plus some arbitrary constant.
If you want to use this indefinite integral / antiderivative to find the definite integral, you need to evaluate it twice: once at the start of your bounds and once at the ends. More information.
One more note is that I believe you swapped the variables analytic and numeric, as the analytic solution is evaluating the function at many points and summing the area, and the numeric solution is based on the integral. Generally people would call the first thing numeric integration and the second thing analytic integration. I've fixed that below.
from scipy.integrate import trapezoid
from scipy.interpolate import CubicSpline
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
def u_time(t_0):
prefactor = 1 / (np.pi**0.25)
def func(time):
return prefactor * np.exp(-(time - t_0) ** 2 / 2)
return func
def u_integral(t_0, times):
u_values = u_time(t_0)(times)
u_abs_2 = np.abs(u_values) ** 2
u_abs_2_int = np.zeros_like(times, dtype=npplex128)
for k in range(1, len(times)):
u_abs_2_int[k] = trapezoid(u_abs_2[:k], times[:k])
return CubicSpline(times, u_abs_2_int)
def gaussian_squared_integral(t_0):
return lambda time: 0.5 * (erf(time - t_0) + erf(t_0))
t_0 = 0
a = -20
b = 100
time_arr = np.linspace(a, b, 1000)
numeric = u_integral(t_0, time_arr)(time_arr)
F_of_x_at_start_point = gaussian_squared_integral(t_0)(a)
analytic = gaussian_squared_integral(t_0)(time_arr) - F_of_x_at_start_point
plt.plot(time_arr, analytic, label="Analytic")
plt.plot(time_arr, numeric, label="Numeric")
plt.legend()
plt.show()
Lastly, I will note that in this loop:
for k in range(1, len(times)):
u_abs_2_int[k] = trapezoid(u_abs_2[:k], times[:k])
If you need to know the value of the integral at each point along the trapezoidal integration, SciPy has a faster alternative here called cumulative_trapezoid, which would let you avoid this loop.