Neka je $f$ po djelovima glatka periodična funkcija s periodom $T=x_1-x_0$ koja na intervalu $[x_0,x_1]$ zadovoljava Dirichletove uvjete:
Tada Fourierov red $$ S(x)=\frac{a_0}{2}+\sum_{n=1}^\infty a_n \cos \frac{2 n \pi}{T}x + b_n \sin \frac{2 n \pi}{T}x, $$ gdje je $$ a_0=\frac{2}{T} \int_{x_0}^{x_1} f(x) \, dx, \\ a_n=\frac{2}{T} \int_{x_0}^{x_1} f(x) \cos \frac{2 n \pi}{T}x \,dx,\\ b_n=\frac{2}{T} \int_{x_0}^{x_1} f(x) \sin \frac{2 n \pi}{T}x \,dx, $$ konvergira u svakoj točki $x\in [x_0,x_1]$ i vrijedi $$ S(x)=\begin{cases} f(x), \textrm{ ako je } f \textrm{ neprekidna u točki } x,\\ \frac{1}{2} [f(x-0)+f(x+0)], \textrm{ inače.} \end{cases} $$
In [1]:
using SymPy
using Plots
Definirajmo simboličku varijablu $x$, simboličku cjelobrojnu (integer) varijablu $n$, i granice intervala $x_0$ i $x_1$.
In [2]:
x=symbols("x",real=true)
n=symbols("n",integer=true)
x₀=symbols("x₀")
x₁=symbols("x₁")
T=symbols("T")
Out[2]:
In [7]:
a(n)=2*integrate(f(x)*cos(2*PI*n*x/T),(x,x₀,x₁))/T
b(n)=2*integrate(f(x)*sin(2*PI*n*x/T),(x,x₀,x₁))/T
Out[7]:
In [5]:
PI
Out[5]:
In [25]:
f(x)=x
x₀=-PI
x₁=PI
T=x₁-x₀
Out[25]:
In [10]:
a(1)
Out[10]:
In [11]:
[a(n) for n=0:5]
Out[11]:
In [12]:
[b(n) for n=0:5]
Out[12]:
In [15]:
x0=Float64(x₀)
x1=Float64(x₁)
K=20
S=a(0)/2+sum([a(n)*cos(2*PI*n*x/T)+b(n)*sin(2*PI*n*x/T) for n=1:K])
g(x)=S(x)
title="Funkcija i "*string(K)*" članova reda"
plot(f,x0,x1)
plot!(g,x0,x1,title=title)
Out[15]:
In [16]:
S
Out[16]:
In [17]:
G₁=series(sin(x),x,0,7)
Out[17]:
In [18]:
G₂=series(cos(x),x,0,7)
Out[18]:
In [19]:
G₃=simplify(G₁*G₂)
Out[19]:
In [20]:
G₃.removeO()
Out[20]:
Sada razvoj u Fourier-ov red
In [21]:
# Ovo učitava funkcije za Fourierovu analizu i Laplaceovu transformaciju
import_from(sympy)
In [26]:
c=fourier_series(x^2,(x,x₀,x₁))
Out[26]:
In [27]:
# Prvih 10 harmonika
c₁=c.truncate(10)
Out[27]:
In [28]:
# Numerička vrijednost u točki x=1
Float64(c₁(1))
Out[28]:
In [29]:
X = range(x0,x1,length=100)
css=[Float64(c₁(x)) for x in X]
Out[29]:
In [30]:
plot(x->x^2,x0,x1)
plot!(X,css)
Out[30]:
In [31]:
f₁=x
K=20
cs=fourier_series(f₁,(x,x₀,x₁)).truncate(K)
css=[Float64(cs(x)) for x in X]
title="Funkcija i "*string(K)*" članova reda"
plot(f,x0,x1)
plot!(X,css,title=title)
Out[31]:
Fourierove koeficijente čemo izračunati numeričkom integracijom (vidi Numeričko integriranje) koristeći Julia naredbu quadgk.
Ovaj dio je izrađen prema bilježnici lecture-2.ipynb Stevena Johnsona izrađenoj za predmet 18.303.
Definirajmo funkcije sinecoef i coscoef koje numerički računaju koeficijente. Parametar abstol je toleranca numeričke integracije: želimo da je greška mala u odnosu na $\sqrt{\int_{x_0}^{x_1} |f(x)|^2 dx}$.
In [32]:
using QuadGK
In [34]:
?quadgk
Out[34]:
In [33]:
Tf=Float64(T)
sinecoef(f, m, x0, x1, Tf) = 2 * quadgk(x -> f(x) * sin(2*m*π*x/Tf)/Tf, x0,x1)[1]
coscoef(f, m, x0, x1, Tf) = 2 * quadgk(x -> f(x) * cos(2*m*π*x/Tf)/Tf, x0,x1)[1]
# i druga funkcija koja računa na vektoru prirodnih brojeva
sinecoef(f, M::AbstractVector,x0,x1,Tf) = Float64[sinecoef(f,m,x0,x1,Tf) for m in M]
coscoef(f, M::AbstractVector,x0,x1,Tf) = Float64[coscoef(f,m,x0,x1,Tf) for m in M]
Out[33]:
In [35]:
# First, define a function to evaluate N terms of the sine series, given the coefficients b
function fouriersum(a, b, x, T)
f = a[1]/2
for n = 1:length(b)
f += a[n+1]* cos(2*n*π*x/T) + b[n] * sin(2*n*π*x/T)
end
return f
end
fouriersum(a, b, X::AbstractVector, T) = Float64[fouriersum(a, b,x,T) for x in X]
Out[35]:
In [36]:
# Uzmimo više točaka (inače "nestane" Gibbsov efekt!!)
X=range(x0,x1,length=1000)
# try for n=1:2:99
n=50
an = coscoef(f, 0:n, x0, x1, Tf)
bn = sinecoef(f, 1:n, x0, x1, Tf)
xlabel="x"
title="Funkcija i "*string(n)*" članova reda"
plot(f,x0,x1,lab="exact f")
plot!(X, fouriersum(an, bn, X, Tf),lab="$n-term sine series",
xlabel=xlabel,title=title)
Out[36]:
In [37]:
f(x)=x^2
x₀=-1
x₁=1
T=x₁-x₀
Out[37]:
In [38]:
a(3)
Out[38]:
In [39]:
f(x)
Out[39]:
In [40]:
a(n)
Out[40]:
In [41]:
b(2)
Out[41]:
In [42]:
f(2)
Out[42]:
In [43]:
plot(f,x₀,x₁)
Out[43]:
In [44]:
x0=Float64(x₀)
x1=Float64(x₁)
# Try for N=1:10
N=10
S=a(0)/2+sum([a(n)*cos(2*PI*n*x/T)+b(n)*sin(2*PI*n*x/T) for n=1:N])
g(x)=S(x)
title="Funkcija i "*string(N)*" članova reda"
plot(f(x),x₀,x₁)
plot!(g,x₀,x₁,title=title)
Out[44]:
In [45]:
S
Out[45]:
In [47]:
# Po dijelovima definirana funkcija
# f(x)=x<0 ? 0 : x
# f(x)=sympy.Piecewise((0,x<=0),(x,x>0))
f(x)=x*Heaviside(x)
plot(f(x),x₀,x₁)
Out[47]:
In [48]:
# Try for N=1:10
N=10
S=a(0)/2+sum([a(n)*cos(2*PI*n*x/T)+b(n)*sin(2*PI*n*x/T) for n=1:N])
g(x)=S(x)
title="Funkcija i "*string(N)*" članova reda"
plot(f(x),x₀,x₁)
plot!(g,x₀,x₁,title=title)
Out[48]:
In [36]:
S
Out[36]: