質問編集履歴
1
コード変更
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -7,52 +7,6 @@
|
|
|
7
7
|
from pylab import *
|
|
8
8
|
|
|
9
9
|
from scipy.integrate import quad
|
|
10
|
-
|
|
11
|
-
|
|
12
|
-
|
|
13
|
-
def fourier(fun, n_max):
|
|
14
|
-
|
|
15
|
-
a = []
|
|
16
|
-
|
|
17
|
-
b = []
|
|
18
|
-
|
|
19
|
-
for n in range(0, 1+n_max):
|
|
20
|
-
|
|
21
|
-
res, err = quad(lambda x:fun(x)*cos(n*x), -pi, pi)
|
|
22
|
-
|
|
23
|
-
a.append(res/pi)
|
|
24
|
-
|
|
25
|
-
for n in range(1, 2+n_max):
|
|
26
|
-
|
|
27
|
-
res, err = quad(lambda x:fun(x)*sin(n*x), pi, pi)
|
|
28
|
-
|
|
29
|
-
b.append(res/pi)
|
|
30
|
-
|
|
31
|
-
def fn(x):
|
|
32
|
-
|
|
33
|
-
sum = a[0] / 2
|
|
34
|
-
|
|
35
|
-
for n in range(1, 1+n_max):
|
|
36
|
-
|
|
37
|
-
sum += a[n]*cos(n*x) + b[n]*sin(n*x)
|
|
38
|
-
|
|
39
|
-
return sum
|
|
40
|
-
|
|
41
|
-
return fn
|
|
42
|
-
|
|
43
|
-
|
|
44
|
-
|
|
45
|
-
def f(x):
|
|
46
|
-
|
|
47
|
-
x = (x + pi) % (pi * 2) - pi
|
|
48
|
-
|
|
49
|
-
if x >= 0:
|
|
50
|
-
|
|
51
|
-
return 1
|
|
52
|
-
|
|
53
|
-
else:
|
|
54
|
-
|
|
55
|
-
return 0
|
|
56
10
|
|
|
57
11
|
|
|
58
12
|
|
|
@@ -66,32 +20,106 @@
|
|
|
66
20
|
|
|
67
21
|
axis([x_min, x_max, y_min, y_max])
|
|
68
22
|
|
|
23
|
+
xs = linspace(x_min, x_max, 256)
|
|
69
24
|
|
|
70
25
|
|
|
71
|
-
f_fn = fourier(f, 10)
|
|
72
26
|
|
|
73
|
-
|
|
27
|
+
n_max = 5
|
|
74
28
|
|
|
75
|
-
plot(xs, amap(f, xs), 'b:', lw=1)
|
|
76
29
|
|
|
30
|
+
|
|
31
|
+
def fun(x):
|
|
32
|
+
|
|
33
|
+
x = (x + pi) % (pi * 2) - pi
|
|
34
|
+
|
|
35
|
+
f = abs(x)
|
|
36
|
+
|
|
37
|
+
|
|
38
|
+
|
|
39
|
+
return f
|
|
40
|
+
|
|
41
|
+
"""
|
|
42
|
+
|
|
43
|
+
if x >= 0:
|
|
44
|
+
|
|
45
|
+
return 1
|
|
46
|
+
|
|
47
|
+
else:
|
|
48
|
+
|
|
49
|
+
return 0
|
|
50
|
+
|
|
51
|
+
"""
|
|
52
|
+
|
|
53
|
+
|
|
54
|
+
|
|
55
|
+
def fourier(fun, n_max):
|
|
56
|
+
|
|
57
|
+
a = []
|
|
58
|
+
|
|
59
|
+
b = []
|
|
60
|
+
|
|
61
|
+
res, err = quad(lambda x:fun(x), -pi, pi)
|
|
62
|
+
|
|
63
|
+
a += [res/pi]
|
|
64
|
+
|
|
65
|
+
b += [0.0]
|
|
66
|
+
|
|
67
|
+
for n in range(1, 1+n_max):
|
|
68
|
+
|
|
69
|
+
res, err = quad(lambda x:fun(x)*cos(n*x), -pi, pi)
|
|
70
|
+
|
|
71
|
+
a+=[res/pi]
|
|
72
|
+
|
|
73
|
+
for n in range(1, 1+n_max):
|
|
74
|
+
|
|
75
|
+
res, err = quad(lambda x:fun(x)*sin(n*x), pi, pi)
|
|
76
|
+
|
|
77
|
+
b+=[res/pi]
|
|
78
|
+
|
|
79
|
+
|
|
80
|
+
|
|
81
|
+
return a, b
|
|
82
|
+
|
|
83
|
+
|
|
84
|
+
|
|
85
|
+
a, b = fourier(fun, n_max)
|
|
86
|
+
|
|
87
|
+
|
|
88
|
+
|
|
89
|
+
def fourier_function(x,a,b):
|
|
90
|
+
|
|
91
|
+
f = 0.5 * a[0]
|
|
92
|
+
|
|
93
|
+
for n in range(1, len(a)):
|
|
94
|
+
|
|
95
|
+
f += a[n]*cos(n*x) + b[n]*sin(n*x)
|
|
96
|
+
|
|
97
|
+
return f
|
|
98
|
+
|
|
99
|
+
|
|
100
|
+
|
|
101
|
+
f_fn = amap(lambda x:fourier_function(x,a,b),xs)
|
|
102
|
+
|
|
103
|
+
|
|
104
|
+
|
|
105
|
+
|
|
106
|
+
|
|
77
|
-
plot(xs, amap(f
|
|
107
|
+
#plot(xs, amap(f, xs), 'b:', lw=1)
|
|
108
|
+
|
|
109
|
+
plot(xs, f_fn, 'r-', lw=1)
|
|
78
110
|
|
|
79
111
|
show()
|
|
80
112
|
|
|
81
113
|
```
|
|
82
114
|
|
|
83
|
-
|
|
84
|
-
|
|
85
|
-
|
|
86
|
-
|
|
87
|
-
理想は以上のようなグラフなのですが、
|
|
88
|
-
|
|
89
|
-
|
|
90
|
-
|
|
91
|
-
実際は以上のようなグラフになってしまいます。
|
|
92
116
|
|
|
93
117
|
|
|
94
118
|
|
|
95
|
-
|
|
119
|
+
a, bをちゃんと求めるようにして実行できるように組んでみたのですが、やはりうまくいきません。
|
|
96
120
|
|
|
121
|
+
|
|
122
|
+
|
|
123
|
+
何度も申し訳ないのですがアドバイスいただけたらありがたいです。
|
|
124
|
+
|
|
97
|
-
よろしくお
|
|
125
|
+
よろしくお願いします。
|