A Fourier series approximation of a discontinuous eddy function

[Remusco]

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:

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:

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)

 

[thomsj4]

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$.

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}")