Fixed gitignore
Added fourier.py program for lab1 More work done on lab1
This commit is contained in:
parent
ba1a333005
commit
90d1f589b5
4 changed files with 208 additions and 23 deletions
2
.gitignore
vendored
2
.gitignore
vendored
|
|
@ -1 +1 @@
|
|||
.pdf
|
||||
*.pdf
|
||||
|
|
|
|||
|
|
@ -34,6 +34,15 @@
|
|||
}
|
||||
})
|
||||
|
||||
#show raw.where(block: true): code => {
|
||||
show raw.line: line => {
|
||||
text(fill: gray)[#line.number]
|
||||
h(1em)
|
||||
line.body
|
||||
}
|
||||
code
|
||||
}
|
||||
|
||||
#show heading: it => {
|
||||
set text(size: 14pt, weight: "bold")
|
||||
block(above: 1.5em, below: 1.5em, it)
|
||||
|
|
@ -60,3 +69,5 @@
|
|||
#doc
|
||||
]
|
||||
|
||||
#let numbered_eq(content) = math.equation(block: true, numbering: "(1)", content)
|
||||
|
||||
|
|
|
|||
27
chastotnie-methods/lab1/fourier.py
Normal file
27
chastotnie-methods/lab1/fourier.py
Normal file
|
|
@ -0,0 +1,27 @@
|
|||
import sys
|
||||
import math
|
||||
|
||||
function_formulas = [
|
||||
# 1. Square wave
|
||||
{
|
||||
'an': lambda n: 3 if n == 0 else 0,
|
||||
'bn': lambda n: 0 if n % 2 == 0 else 2 / (math.pi * n),
|
||||
'cn': lambda n: 3 if n == 0 else (0 if n % 2 == 0 else 1j / (math.pi * n))
|
||||
}
|
||||
]
|
||||
|
||||
index, nn = int(sys.argv[1]), int(sys.argv[2])
|
||||
|
||||
formulas = function_formulas[index - 1]
|
||||
coefficients = {'an': [], 'bn': [], 'cn': []}
|
||||
|
||||
for coefficient in coefficients:
|
||||
for n in range(0, nn + 1):
|
||||
calc = formulas[coefficient](n)
|
||||
coefficients[coefficient].append(f'{n}: {calc:.5g}')
|
||||
if coefficient == 'cn' and n != 0:
|
||||
calc = formulas[coefficient](-n)
|
||||
coefficients[coefficient].insert(0, f'{-n}: {calc:.5g}')
|
||||
|
||||
for coefficient in coefficients:
|
||||
print(f'{coefficient} values: ', coefficients[coefficient])
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
#import "../lab-template.typ"
|
||||
#import lab-template: lab
|
||||
#import lab-template: *
|
||||
#import "@preview/cetz:0.2.2": canvas, plot
|
||||
|
||||
#show: doc => lab(
|
||||
|
|
@ -9,7 +9,7 @@
|
|||
)
|
||||
|
||||
= Вещественные функции
|
||||
Придумайте числа $a, b, t_0, t_1, t_2$ такие, что $a, b > 0$ и $t_2 > t_1 > t_0$. Рассмотрите следующие функции $f : RR -> RR$ и для каждой из функции:
|
||||
Придумайте числа $a, b, t_0, t_1, t_2$ такие, что $a, b > 0$ и $t_2 > t_1 > t_0 > 0$. Рассмотрите следующие функции $f : RR -> RR$ и для каждой из функции:
|
||||
|
||||
- Постройте график $f(t)$
|
||||
|
||||
|
|
@ -32,26 +32,26 @@ $ f(t) = cases(
|
|||
b\, & t in [t_1, t_2)
|
||||
) $
|
||||
|
||||
Пусть $a = 1, b = 2$ и $t_0 = 0, t_1 = 1, t_2 = 2$. В таком случае конечная функция будет вида:
|
||||
Пусть $a = 1, b = 2$ и $t_0 = 1, t_1 = 2, t_2 = 3$. В таком случае конечная функция будет вида:
|
||||
|
||||
$ f(t) = cases(
|
||||
1\, & t in [0, 1),
|
||||
2\, & t in [1, 2)
|
||||
1\, & t in [1, 2),
|
||||
2\, & t in [2, 3)
|
||||
) $
|
||||
|
||||
Построим график данной функции:
|
||||
|
||||
#figure(
|
||||
canvas(length: 1.75cm, {
|
||||
plot.plot(size: (3, 2),
|
||||
canvas(length: 1.25cm, {
|
||||
plot.plot(size: (8, 4),
|
||||
x-tick-step: 1, y-tick-step: 1,
|
||||
x-min: -3, x-max: 3,
|
||||
y-min: -1, y-max: 4,
|
||||
y-min: 0, y-max: 3,
|
||||
x-grid: true, y-grid: true,
|
||||
x-label: [$t$], y-label: [$f(t)$],
|
||||
{
|
||||
plot.add(domain: (-2.9, 2.9), t => {
|
||||
let (t0, t1, t2) = (0, 1, 2)
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
let (t0, t1, t2) = (1, 2, 3)
|
||||
if t < 0 { t = -t + 1 }
|
||||
let abs_rem = calc.rem(t, t2 - t0)
|
||||
if abs_rem >= t0 and abs_rem < t1 {
|
||||
|
|
@ -59,27 +59,174 @@ $ f(t) = cases(
|
|||
} else {
|
||||
return 2
|
||||
}
|
||||
}, line: "hvh")
|
||||
}, line: "hvh", samples: 100)
|
||||
}
|
||||
)
|
||||
}), caption: [График квадратной волны $f(t)$]
|
||||
) <sq_wave>
|
||||
#pagebreak()
|
||||
|
||||
Период функции равен $T = t_2 - t_0 = 2$, рассмотрим функцию на промежутке $[-T/2; T/2]$ и найдём коэффициенты $a_0, a_n, b_n, c_n$:
|
||||
Период функции равен $T = t_2 - t_0 = 2$, рассмотрим функцию на промежутке $[1; 3]$ и найдём коэффициенты $a_0, a_n, b_n, c_n$:
|
||||
|
||||
$ a_0 &= 1/T integral^(T/2)_(-T/2) f(t)upright(d)t = 1/2 integral^1_(-1) f(t)upright(d)t = 1/2 (integral^0_(-1) 2 + integral^1_0 1)upright(d)t \
|
||||
&= 1/2 (2t bar.v^0_(-1) + t bar.v^1_0) = 1/2 (2 + 1) = 3/2 $
|
||||
$ a_0 &= 1/(T/2) integral^(3)_(1) f(t)upright(d)t = integral^3_1 f(t)upright(d)t = (integral^2_1 1 + integral^3_2 2)upright(d)t = (lr(t|)^2_1 + lr(2t|)^3_2) \
|
||||
&= 2 - 1 + 6 - 4 = 3 $
|
||||
|
||||
$ a_n &= 2/T integral^(T/2)_(-T/2) f(t)cos(omega_n t)upright(d)t = integral^(1)_(-1) f(t)cos(pi n t)upright(d)t \
|
||||
&= (integral^0_(-1) 2cos(pi n t) + integral^(1)_0 cos(pi n t))upright(d)t \
|
||||
&= (2 sin(pi n t))/(pi n) bar.v^0_(-1) + sin(pi n t)/(pi n) bar.v^(1)_0 \
|
||||
&= -(2 sin(-pi n))/(pi n) + sin(pi n)/(pi n) = (3 sin(pi n))/(pi n) = 0 $
|
||||
$ a_n &= 1/(T/2) integral^3_1 f(t)cos(omega_n t)upright(d)t = integral^3_1 f(t)cos(pi n t)upright(d)t \
|
||||
&= (integral^2_1 cos(pi n t) + integral^3_2 2cos(pi n t))upright(d)t = lr(sin(pi n t)/(pi n)|)^2_1 + lr((2sin(pi n t))/(pi n)|)^3_2 \
|
||||
&= sin(2pi n)/(pi n) - sin(pi n)/(pi n) + (2sin(3pi n))/(pi n) - (2sin(2pi n))/(pi n) \
|
||||
&= 0 "(т.к. " sin(3pi n) = sin(2pi n) = sin(pi n) = 0")" $
|
||||
|
||||
$ b_n &= 2/T integral^(T/2)_(-T/2) f(t)sin(omega_n t)upright(d)t = integral^(1)_(-1) f(t)sin(pi n t)upright(d)t \
|
||||
&= (integral^0_(-1) 2sin(pi n t) + integral^(1)_0 sin(pi n t))upright(d)t \
|
||||
&= -(2 cos(pi n t))/(pi n) bar.v^0_(-1) - cos(pi n t)/(pi n) bar.v^(1)_0 \
|
||||
&= (2 cos(-pi n))/(pi n) - cos(pi n)/(pi n) = cos(pi n)/(pi n) $
|
||||
$ b_n &= 1/(T/2) integral^3_1 f(t)sin(omega_n t)upright(d)t = integral^3_1 f(t)sin(pi n t)upright(d)t \
|
||||
&= (integral^2_1 sin(pi n t) + integral^3_2 2sin(pi n t))upright(d)t = lr(-cos(pi n t)/(pi n)|)^2_1 - lr((2cos(pi n t))/(pi n)|)^3_2 \
|
||||
&= -cos(2pi n)/(pi n) + cos(pi n)/(pi n) - (2cos(3pi n))/(pi n) + (2cos(2pi n))/(pi n) \
|
||||
&= cos(2pi n)/(pi n) + cos(pi n)/(pi n) - (2cos(pi n))/(pi n) = 1/(pi n) - cos(pi n)/(pi n) = (1 - (-1)^n)/(pi n) \
|
||||
&= cases(
|
||||
0\, & n = 2k,
|
||||
2/(pi n)\, & n = 2k - 1
|
||||
), k in ZZ $
|
||||
|
||||
$ c_n &= 1/T integral^3_1 f(t)e^(i omega_n t) = 1/2 integral^3_1 f(t)e^(i pi n t) = 1/2 (integral^2_1 e^(i pi n t) + integral^3_2 2e^(i pi n t))upright(d)t \
|
||||
&= 1/2 (lr(-(i e^(i pi n t))/(pi n)|)^2_1 - lr((2i e^(i pi n t))/(pi n)|)^3_2) = -(i e^(2i pi n))/(2pi n) + (i e^(i pi n))/(2pi n) - (2i e^(3i pi n))/(2pi n) + (2i e^(2i pi n))/(2pi n) \
|
||||
&= (i e^(2i pi n))/(2pi n) - (i e^(i pi n))/(2pi n) = i/(2pi n) - (i(-1)^n)/(2pi n) = (i(1 - (-1)^n))/(2pi n) \
|
||||
&= cases(
|
||||
0\, & n = 2k,
|
||||
i/(pi n)\, & n = 2k - 1
|
||||
), k in ZZ $
|
||||
|
||||
$ a_0 &= c_0 = 3, a_n = 0, \
|
||||
b_n &= cases(
|
||||
0\, & n = 2k,
|
||||
2/(pi n)\, & n = 2k - 1
|
||||
), k in ZZ \
|
||||
c_n &= cases(
|
||||
0\, & n = 2k,
|
||||
i/(pi n)\, & n = 2k - 1
|
||||
), k in ZZ $
|
||||
|
||||
Вычислим значения коэффициентов для $n = 0, 1, 2$:
|
||||
|
||||
$ A = vec(a_0, a_1, a_2) = vec(3, 0, 0), B = vec(b_0, b_1, b_2) = vec(0, 2/pi, 0), C = vec(c_0, c_1, c_2) = vec(3, i/pi, 0) $
|
||||
|
||||
Вручную считать тяжко, поэтому напишем программу на _Python_, которая будет принимать на вход $N$ и считать коэффициенты по найденным формулам:
|
||||
|
||||
#figure(```python
|
||||
import sys
|
||||
import math
|
||||
|
||||
function_formulas = [
|
||||
# 1. Square wave
|
||||
{
|
||||
'an': lambda n: 3 if n == 0 else 0,
|
||||
'bn': lambda n: 0 if n % 2 == 0 else 2 / (math.pi * n),
|
||||
'cn': lambda n: 3 if n == 0 else (0 if n % 2 == 0 else 1j / (math.pi * n))
|
||||
|
||||
}
|
||||
]
|
||||
|
||||
index, nn = int(sys.argv[1]), int(sys.argv[2])
|
||||
|
||||
formulas = function_formulas[index - 1]
|
||||
coefficients = {'an': [], 'bn': [], 'cn': []}
|
||||
|
||||
for coefficient in coefficients:
|
||||
for n in range(0, nn + 1):
|
||||
calc = formulas[coefficient](n)
|
||||
coefficients[coefficient].append(f'{n}: {calc:.5g}')
|
||||
|
||||
for coefficient in coefficients:
|
||||
print(f'{coefficient} values: ', coefficients[coefficient])
|
||||
```, caption: "Код программы для подсчета коэффициентов")
|
||||
|
||||
Программа принимает на вход номер задания соответствующей функции и произвольное $N$ и выводит подсчитанные коэффициенты $a_n, b_n, c_n$ для $n = 0,1,2...N$. Список _function_formulas_ хранит формулы для расчетов коэффитциентов в виде лямда выражений для каждой функции, данный список будет расширяться по мере нахождения формул для других функций.
|
||||
|
||||
Воспользуемся программой для нахождения коэффициентов для $N = 3$, нас интересует функция 1 - квадратная волна, вызывыаем программу с аргументами _1_ и _3_:
|
||||
|
||||
#figure(```
|
||||
$ python fourier.py 1 3
|
||||
an values: ['0: 3', '1: 0', '2: 0', '3: 0']
|
||||
bn values: ['0: 0', '1: 0.63662', '2: 0', '3: 0.21221']
|
||||
cn values: ['-3: -0-0.1061j', '-2: 0', '-1: -0-0.31831j', '0: 3', '1: 0+0.31831j', '2: 0', '3: 0+0.1061j']
|
||||
```, caption: "Ввывод программы для первой функции") <prog>
|
||||
|
||||
Значения коэффициентов для $n = 0, 1, 2$ совпадают с посчитанными вручную, значит высока вероятность, что программа работает корректно. Также стоит заметить, что программа округляет коэффициенты до 5 знаков после запятой и считает дополнительно коэффициенты $c_n$ для отрицательных $n$
|
||||
|
||||
Запишем частичные суммы Фурье $F_3$ и $G_3$ с учётом подсчитанных коэффициентов:
|
||||
|
||||
$ F_3(t) = 3/2 + 0.63662sin(pi t) + 0.21221sin(3pi t) $
|
||||
$ G_3(t) = -0.1061i e^(-3i pi t) -0.31831i e^(-i pi t) + 3/2 + 0.31831i e^(i pi t) + 0.1061i e^(3i pi t) $
|
||||
|
||||
С помощью программы из @prog, подсчитаем также коэффициенты для $N = 5, 7, 9$, после чего построим график $F_N (t)$ для значений $N = 1, 3, 5, 7, 9$ и сравним их с $f(t)$. Графики $G_N (t)$ получатся аналогичны $F_N (t)$, поскольку изначальная функция не является комплексной (вещественная часть будет совпадать, а комплексной части - нет):
|
||||
|
||||
#figure(
|
||||
canvas(length: 1.75cm, {
|
||||
plot.plot(size: (8, 4),
|
||||
x-tick-step: 1, y-tick-step: 1,
|
||||
x-min: -3, x-max: 3,
|
||||
y-min: 0, y-max: 3,
|
||||
x-grid: true, y-grid: true,
|
||||
x-label: [$t$], y-label: [$F_N (t)$],
|
||||
{
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
let (t0, t1, t2) = (1, 2, 3)
|
||||
if t < 0 { t = -t + 1 }
|
||||
let abs_rem = calc.rem(t, t2 - t0)
|
||||
if abs_rem >= t0 and abs_rem < t1 {
|
||||
return 1
|
||||
} else {
|
||||
return 2
|
||||
}
|
||||
}, line: "hvh", samples: 100, label: [$f(t)$#v(1em)])
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
3/2 + 0.63662*calc.sin(calc.pi*t)
|
||||
}, line: "spline", samples: 100, label: [$F_1(t)$#v(1em)])
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
3/2 + 0.63662*calc.sin(calc.pi*t) + 0.21221*calc.sin(3*calc.pi*t)
|
||||
}, line: "spline", samples: 100, label: [$F_3(t)$#v(1em)])
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
3/2 + 0.63662*calc.sin(calc.pi*t) + 0.21221*calc.sin(3*calc.pi*t) + 0.12732*calc.sin(5*calc.pi*t)
|
||||
}, line: "spline", samples: 100, label: [$F_5(t)$#v(1em)])
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
3/2 + 0.63662*calc.sin(calc.pi*t) + 0.21221*calc.sin(3*calc.pi*t) + 0.12732*calc.sin(5*calc.pi*t) + 0.090946*calc.sin(7*calc.pi*t)
|
||||
}, line: "spline", samples: 100, label: [$F_7(t)$#v(1em)])
|
||||
plot.add(domain: (-2.99, 2.99), t => {
|
||||
3/2 + 0.63662*calc.sin(calc.pi*t) + 0.21221*calc.sin(3*calc.pi*t) + 0.12732*calc.sin(5*calc.pi*t) + 0.090946*calc.sin(7*calc.pi*t) + 0.070736*calc.sin(9*calc.pi*t)
|
||||
}, line: "spline", samples: 100, style: (stroke: purple), label: [$F_9(t)$])
|
||||
}
|
||||
)
|
||||
}), caption: [Графики частичных сумм Фурье $F_N (t)$ и $f(t)$]
|
||||
) <sq_wave_F>
|
||||
|
||||
Из графиков видно, что сумма Фурье $F_N (t)$ сходится к исходной функции $f(t)$ при $N -> infinity$ - коэффициенты $a_n, b_n, c_n$ подсчитаны верно.
|
||||
|
||||
Проверим выполнение равенства Парсеваля:
|
||||
|
||||
$ integral^3_1 (f(t))^2 upright(d)t = a_0^2/2 + sum^infinity_(n=1) (a_n^2 + b_n^2) $
|
||||
|
||||
Вычислим левую часть:
|
||||
|
||||
$ integral^3_1 (f(t))^2 upright(d)t = (integral^2_1 1^2 + integral^3_2 2^2) upright(d)t = lr(t|)^2_1 + lr(4t|)^3_2 = 2 - 1 + 12 - 8 = 5 $
|
||||
|
||||
#pagebreak()
|
||||
|
||||
Подробнее рассмотрим ряд из правой части равенства:
|
||||
|
||||
$ sum^infinity_(n=1) (a_n^2 + b_n^2) = sum^infinity_(n=1) b_n^2 $
|
||||
|
||||
Коэффициент $b_n$ задан частями, при чётных $n$ он равен $0$, а при нечётных $2/(pi n)$. В таком случае можно игнорировать все чётные значения $n$ ($2k$) и рассматривать лишь нечётные ($2k - 1$):
|
||||
|
||||
$ sum^infinity_(n=1) b_n^2 = sum^infinity_(k=1) (2/(pi(2k - 1)))^2 = 4/pi^2 sum^infinity_(k=1) 1/(2k - 1)^2 $
|
||||
|
||||
Полученный нами ряд состоит из элементов ${1, 1/9, 1/25, 1/49...}$. Рассмотрим ряды ${1/k^2; k in NN}$ ${1/(4k^2); k in NN}$. Они состоят соответсвенно из элементов ${1, 1/4, 1/9, 1/16, 1/25, 1/36, 1/49...}$ и ${1/4, 1/16, 1/36...}$. Заметим, что исключив из ряда ${1/k^2; k in NN}$ элементы ряда ${1/(4k^2); k in NN}$, мы получим наш исходный ряд, следовательно, сумма исходного ряда будет равна разности сумм введённых рядов:
|
||||
|
||||
$ sum^infinity_(k=1) 1/(2k - 1)^2 = sum^infinity_(k=1) 1/k^2 - sum^infinity_(k=1) 1/(4k^2) = 3/4 sum^infinity_(k=1) 1/k^2 = 3/4 dot pi^2/6 = pi^2/8 $
|
||||
|
||||
$ sum^infinity_(n=1) (a_n^2 + b_n^2) = sum^infinity_(n=1) b_n^2 = 4/pi^2 dot pi^2/8 = 1/2 $
|
||||
|
||||
$ a_0^2/2 + sum^infinity_(n=1) (a_n^2 + b_n^2) = 9/2 + 1/2 = 10/2 = 5 $
|
||||
|
||||
$ underbrace(integral^3_1 (f(t))^2 upright(d)t, 5) = underbrace(a_0^2/2 + sum^infinity_(n=1) (a_n^2 + b_n^2), 5) $
|
||||
|
||||
Левая и правая часть равны - равенство Парсеваля выполняется.
|
||||
|
||||
2. Любая *чётная* периодическая функция по вашему выбору.
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue