A Fourier series approximation of a discontinuous eddy function

AI Thread Summary
The discussion focuses on approximating a discontinuous velocity function U_theta(r) using a Fourier cosine series. The function is defined piecewise, with different expressions for r values less than and greater than a specified radius R. The Fourier coefficients are derived through integration over two intervals, leading to expressions for a_0 and a_n, which decay as 1/n^2. The implementation uses SciPy for numerical integration and includes a plot comparing the original function with its Fourier series approximation. The results indicate that increasing the parameter L improves the approximation accuracy, with a maximum error calculated for a specified number of terms in the series.
Remusco
Messages
30
Reaction score
3
TL;DR Summary
I have a velocity function of radius r, that I would like to write a Fourier Series expansion to approximate.
I have a velocity function U_theta(r), where r is radius from the origin. This function is defined as:

def U_theta(r):
s=2.5
R=1*s
U=1
if r<=R:
U_theta=r*U/R
if r>R:
U_theta=R*U/r
return U_theta

The plot of U_theta as a function of R is shown as:
1738207508318.png


How can I write a Fourier series to approximate this curve?
I have a start on this, but am not getting a satisfactory result, as shown below:
1738207970192.png


My code that does this:

from scipy.integrate import quad
import numpy as np

def U_theta(r):
s=2.5
R=1*s
U=1
if r<=R:
U_theta=r*U/R
if r>R:
U_theta=R*U/r
return U_theta

L=5
def fourier(r):
A0=(1/L)*quad(lambda r: U_theta(r),0,L)[0]

An=0
Bn=0
for n in range(100):
An+=np.cos(n*np.pi*r/L)*(2/L)*quad(lambda r: U_theta(r)*np.cos(n*np.pi*r/L),0,L)[0]
Bn+=np.sin(n*np.pi*r/L)*(2/L)*quad(lambda r: U_theta(r)*np.sin(n*np.pi*r/L),0,L)[0]
return (A0/2)+An+Bn

fourier_x_vals=np.linspace(0,L,200)
fourier_y_vals=fourier(fourier_x_vals)
plt.plot(fourier_x_vals,fourier_y_vals)

U_theta_y_vals=[U_theta(r) for r in fourier_x_vals]
plt.plot(fourier_x_vals,U_theta_y_vals)
 
Mathematics news on Phys.org
Because ##U_\theta(r)## can be extended evenly around ##r=0## (##U_\theta(-r) = U_\theta(r)##), a Fourier cosine series is used. This means sine terms are zero (##B_n = 0##).

We work on ##[0, L]##, with ##L## significantly larger than ##R## (e.g., ##L = 4R##, ##5R##, or ##10R##).

The Fourier cosine series is:
$$U_\theta(r) \approx \frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \cos(\frac{n \pi r}{L})$$

The coefficients are:
$$a_0 = \frac{2}{L} \int_0^L U_\theta(r) dr$$
$$a_n = \frac{2}{L} \int_0^L U_\theta(r)\cos(\frac{n \pi r}{L})dr$$

We split the integrals:

* From ##0## to ##R##: ##U_\theta(r) = \frac{r}{R}## (since ##U=1##).

$$I_1 = \int_0^R \frac{r}{R}\cos(\frac{n\pi r}{L})dr$$

* From ##R## to ##L##: ##U_\theta(r) = \frac{R}{r}##.

$$I_2 = \int_R^L \frac{R}{r}\cos(\frac{n\pi r}{L})dr$$

These are solved exactly. ##I_1## uses integration by parts. ##I_2## involves ##Ci(x) = -\int_x^\infty \frac{\cos t}{t} dt##:
$$I_1 = \frac{L^2}{R(n\pi)^2}(\cos(\frac{n\pi R}{L}) - 1) + \frac{L}{n\pi}\sin(\frac{n\pi R}{L})$$
$$I_2 = R[\text{Ci}(n\pi) - \text{Ci}(\frac{n\pi R}{L})]$$

For ##n \ge 1##:
$$a_n = \frac{2}{L} (I_1 + I_2) = \frac{2L}{(n\pi)^2 R}(\cos(\frac{n\pi R}{L}) - 1) + \frac{2}{n\pi}\sin(\frac{n\pi R}{L}) + \frac{2R}{L}[\text{Ci}(n\pi) - \text{Ci}(\frac{n\pi R}{L})]$$

For ##n=0##:
$$a_0 = \frac{2}{L} \int_0^L U_\theta(r)dr = \frac{R}{L}(1 + 2\ln(\frac{L}{R}))$$

, giving a constant term of ##\frac{a_0}{2}##. The ##a_n## values decay as ##1/n^2##.

Code:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sici  # For Cosine integral

# Parameters (choose a larger L)
R = 2.5
L = 10 * R # Increased L for better approximation
U = 1.0

# Original function
def U_theta(r):
    return (U * r / R) if r <= R else (U * R / r)

# Number of terms in Fourier series
N = 50

# Calculate Fourier coefficients
a_cos = np.zeros(N + 1)
a_cos[0] = (R / L) * (1 + 2 * np.log(L / R))  # a_0

for n in range(1, N + 1):
    k = n * np.pi / L
    part1 = (np.cos(k * R) - 1) * L**2 / (R * (n * np.pi)**2)
    part2 = (L / (n * np.pi)) * np.sin(k * R)
    _, Ci_npi = sici(n * np.pi)
    _, Ci_term = sici(n * np.pi * R / L)
    part3 = R * (Ci_npi - Ci_term)
    a_cos[n] = (2.0 / L) * (part1 + part2 + part3)

# Generate points for plotting
r_vals = np.linspace(0, 4*R, 500)  # Plot up to 4*R
U_exact = np.array([U_theta(r) for r in r_vals])

# Calculate the Fourier series approximation
U_series = np.full_like(U_exact, a_cos[0] / 2)  # Start with a_0/2
for n in range(1, N + 1):
    U_series += a_cos[n] * np.cos(n * np.pi * r_vals / L)

# Plotting
plt.figure(figsize=(8, 6))
plt.plot(r_vals, U_exact, label='Original $U_\\theta(r)$', linewidth=2)
plt.plot(r_vals, U_series, '--', label=f'Fourier Series (N={N})', linewidth=2)
plt.axvline(x=R, color='gray', linestyle=':', label='$r=R$')
plt.xlabel('$r$')
plt.ylabel('$U_\\theta(r)$')
plt.title('Fourier Cosine Series Approximation')
plt.legend()
plt.grid(True)
plt.xlim(0, 4*R) # Limit x-axis
plt.show()

#Calculate error
error = np.abs(U_exact - U_series)
max_error = np.max(error)
print(f"Maximum error with N={N}: {max_error:.4f}")
 
Last edited:
Suppose ,instead of the usual x,y coordinate system with an I basis vector along the x -axis and a corresponding j basis vector along the y-axis we instead have a different pair of basis vectors ,call them e and f along their respective axes. I have seen that this is an important subject in maths My question is what physical applications does such a model apply to? I am asking here because I have devoted quite a lot of time in the past to understanding convectors and the dual...
Fermat's Last Theorem has long been one of the most famous mathematical problems, and is now one of the most famous theorems. It simply states that the equation $$ a^n+b^n=c^n $$ has no solutions with positive integers if ##n>2.## It was named after Pierre de Fermat (1607-1665). The problem itself stems from the book Arithmetica by Diophantus of Alexandria. It gained popularity because Fermat noted in his copy "Cubum autem in duos cubos, aut quadratoquadratum in duos quadratoquadratos, et...
Insights auto threads is broken atm, so I'm manually creating these for new Insight articles. In Dirac’s Principles of Quantum Mechanics published in 1930 he introduced a “convenient notation” he referred to as a “delta function” which he treated as a continuum analog to the discrete Kronecker delta. The Kronecker delta is simply the indexed components of the identity operator in matrix algebra Source: https://www.physicsforums.com/insights/what-exactly-is-diracs-delta-function/ by...

Similar threads

Replies
15
Views
1K
Replies
3
Views
2K
Replies
2
Views
2K
Replies
6
Views
2K
Replies
4
Views
3K
Replies
3
Views
2K
Replies
0
Views
1K
Replies
1
Views
1K
Back
Top