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:
Seemingly by some mathematical coincidence, a hexagon of sides 2,2,7,7, 11, and 11 can be inscribed in a circle of radius 7. The other day I saw a math problem on line, which they said came from a Polish Olympiad, where you compute the length x of the 3rd side which is the same as the radius, so that the sides of length 2,x, and 11 are inscribed on the arc of a semi-circle. The law of cosines applied twice gives the answer for x of exactly 7, but the arithmetic is so complex that the...
Thread 'Video on imaginary numbers and some queries'
Hi, I was watching the following video. I found some points confusing. Could you please help me to understand the gaps? Thanks, in advance! Question 1: Around 4:22, the video says the following. So for those mathematicians, negative numbers didn't exist. You could subtract, that is find the difference between two positive quantities, but you couldn't have a negative answer or negative coefficients. Mathematicians were so averse to negative numbers that there was no single quadratic...
Thread 'Unit Circle Double Angle Derivations'
Here I made a terrible mistake of assuming this to be an equilateral triangle and set 2sinx=1 => x=pi/6. Although this did derive the double angle formulas it also led into a terrible mess trying to find all the combinations of sides. I must have been tired and just assumed 6x=180 and 2sinx=1. By that time, I was so mindset that I nearly scolded a person for even saying 90-x. I wonder if this is a case of biased observation that seeks to dis credit me like Jesus of Nazareth since in reality...

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