A Fourier series approximation of a discontinuous eddy function

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:
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...
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...
I'm interested to know whether the equation $$1 = 2 - \frac{1}{2 - \frac{1}{2 - \cdots}}$$ is true or not. It can be shown easily that if the continued fraction converges, it cannot converge to anything else than 1. It seems that if the continued fraction converges, the convergence is very slow. The apparent slowness of the convergence makes it difficult to estimate the presence of true convergence numerically. At the moment I don't know whether this converges or not.

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