回答編集履歴
3
整形
answer
CHANGED
|
@@ -71,87 +71,75 @@
|
|
|
71
71
|
import numpy as np
|
|
72
72
|
import matplotlib.pyplot as plt
|
|
73
73
|
|
|
74
|
+
|
|
74
75
|
# 周期3次スプライン曲線のパラメータを求める。
|
|
75
76
|
def calculatePeriodicSplineParameter(x):
|
|
76
77
|
# t = 0,1,2,...,x.size-1とする。
|
|
77
|
-
n = x.size-1
|
|
78
|
+
n = x.size - 1
|
|
78
79
|
h = np.ones(n)
|
|
79
80
|
a = x
|
|
80
81
|
|
|
81
82
|
# 対角成分
|
|
82
|
-
Adiag = np.
|
|
83
|
+
Adiag = np.zeros((n, n))
|
|
83
|
-
Adiag[1:,1:]
|
|
84
|
+
Adiag[1:, 1:] = np.diag(2 * (h[:-1] + h[1:]))
|
|
84
|
-
Adiag[0][0]
|
|
85
|
+
Adiag[0][0] = 2 * (h[-1] + h[0])
|
|
85
|
-
|
|
86
|
-
# 対角成分の一つ下の斜めの成分
|
|
87
|
-
Alower = np.
|
|
86
|
+
Alower = np.diag(h[:-1], k=-1)
|
|
88
|
-
Alower[1:,:-1]*=h[:-1]
|
|
89
|
-
|
|
90
|
-
# 対角成分の一つ上の斜めの成分
|
|
91
|
-
Aupper = np.
|
|
87
|
+
Aupper = np.diag(h[:-1], k=1)
|
|
92
|
-
|
|
88
|
+
|
|
93
89
|
A = Adiag + Alower + Aupper
|
|
94
|
-
|
|
95
90
|
A[-1][0] = h[-1]
|
|
96
91
|
A[0][-1] = h[-1]
|
|
97
92
|
|
|
98
|
-
#
|
|
93
|
+
# 対角成分
|
|
99
|
-
|
|
94
|
+
Vdiag = np.zeros((n, n))
|
|
95
|
+
Vdiag[1:, 1:] = np.diag(-(1 / h[:-1] + 1 / h[1:]))
|
|
96
|
+
Vdiag[0][0] = -(1 / h[-1] + 1 / h[0])
|
|
97
|
+
Vlower = np.diag(1 / h[:-1], k=-1)
|
|
98
|
+
Vupper = np.diag(1 / h[:-1], k=1)
|
|
100
99
|
|
|
101
|
-
# 対角成分
|
|
102
|
-
Vdiag = np.eye(n)
|
|
103
|
-
Vdiag[1:,1:] *= -(1/h[:-1]+1/h[1:])
|
|
104
|
-
Vdiag[0][0] *= -(1/h[-1]+1/h[0])
|
|
105
|
-
|
|
106
|
-
# 対角成分の一つ下の斜めの成分
|
|
107
|
-
Vlower = np.eye(n, k=-1)
|
|
108
|
-
Vlower[1:,:-1]*=1/h[:-1]
|
|
109
|
-
|
|
110
|
-
# 対角成分の一つ上の斜めの成分
|
|
111
|
-
Vupper = np.eye(n, k=1)
|
|
112
|
-
Vupper[:-1,1:]*=1/h[:-1]
|
|
113
100
|
V = Vdiag + Vlower + Vupper
|
|
114
|
-
|
|
115
101
|
V[-1][0] = h[-1]
|
|
116
102
|
V[0][-1] = h[-1]
|
|
117
|
-
|
|
103
|
+
|
|
118
|
-
r = 3*np.dot(V,a[:-1])
|
|
104
|
+
r = 3 * np.dot(V, a[:-1])
|
|
119
|
-
|
|
105
|
+
|
|
120
|
-
c = np.
|
|
106
|
+
c = np.linalg.solve(A, r)
|
|
121
|
-
|
|
107
|
+
|
|
122
108
|
c = np.insert(c, c.size, c[0])
|
|
109
|
+
|
|
123
|
-
d = (c[1:]-c[:-1])/(3*h)
|
|
110
|
+
d = (c[1:] - c[:-1]) / (3 * h)
|
|
124
|
-
b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
|
|
111
|
+
b = (a[1:] - a[:-1]) / h - h / 3 * (2 * c[:-1] + c[1:])
|
|
125
112
|
a = a[:-1]
|
|
126
113
|
c = c[:-1]
|
|
127
|
-
return [a,b,c,d]
|
|
114
|
+
return [a, b, c, d]
|
|
128
115
|
|
|
116
|
+
|
|
129
117
|
# 三次関数の係数から, 元の点から間の点を求める。
|
|
130
|
-
def calculateInterpolationPoints(splineParameter,
|
|
118
|
+
def calculateInterpolationPoints(splineParameter,
|
|
119
|
+
interpolationPointsNumber=1000):
|
|
131
120
|
a, b, c, d = splineParameter
|
|
132
|
-
assert(a.size==b.size==c.size==d.size)
|
|
121
|
+
assert (a.size == b.size == c.size == d.size)
|
|
133
122
|
index = np.repeat(np.arange(a.size), interpolationPointsNumber, axis=0)
|
|
134
123
|
t = np.tile(np.linspace(0, 1, interpolationPointsNumber), a.size)
|
|
135
|
-
return a[index] + b[index]*t + c[index]*t**2 + d[index]*t**3
|
|
124
|
+
return a[index] + b[index] * t + c[index] * t**2 + d[index] * t**3
|
|
136
125
|
|
|
126
|
+
|
|
137
127
|
# 補間したい点の座標
|
|
138
|
-
x = np.random.rand(
|
|
128
|
+
x = np.random.rand(100)
|
|
139
129
|
y = np.random.rand(x.size)
|
|
140
|
-
x = np.insert(x,x.size,x[0])
|
|
130
|
+
x = np.insert(x, x.size, x[0])
|
|
141
|
-
y = np.insert(y,y.size,y[0])
|
|
131
|
+
y = np.insert(y, y.size, y[0])
|
|
142
132
|
|
|
143
133
|
parameterX = calculatePeriodicSplineParameter(x)
|
|
144
134
|
parameterY = calculatePeriodicSplineParameter(y)
|
|
145
135
|
|
|
146
|
-
|
|
147
136
|
X = calculateInterpolationPoints(parameterX)
|
|
148
137
|
Y = calculateInterpolationPoints(parameterY)
|
|
149
138
|
|
|
150
139
|
plt.plot(X, Y, label="spline curve")
|
|
151
140
|
plt.plot(x, y, marker='o', linestyle='None', label='original points')
|
|
152
141
|
plt.legend()
|
|
153
|
-
plt.savefig("graph.png", dpi=200,figsize=(4,3))
|
|
142
|
+
plt.savefig("graph.png", dpi=200, figsize=(4, 3))
|
|
154
|
-
|
|
155
|
-
|
|
143
|
+
|
|
156
144
|
```
|
|
157
145
|
|
2
コード修正
answer
CHANGED
|
@@ -98,7 +98,20 @@
|
|
|
98
98
|
# Aの逆行列
|
|
99
99
|
Ainv = np.linalg.inv(A)
|
|
100
100
|
|
|
101
|
+
# 対角成分
|
|
102
|
+
Vdiag = np.eye(n)
|
|
103
|
+
Vdiag[1:,1:] *= -(1/h[:-1]+1/h[1:])
|
|
104
|
+
Vdiag[0][0] *= -(1/h[-1]+1/h[0])
|
|
105
|
+
|
|
106
|
+
# 対角成分の一つ下の斜めの成分
|
|
107
|
+
Vlower = np.eye(n, k=-1)
|
|
108
|
+
Vlower[1:,:-1]*=1/h[:-1]
|
|
109
|
+
|
|
110
|
+
# 対角成分の一つ上の斜めの成分
|
|
101
|
-
|
|
111
|
+
Vupper = np.eye(n, k=1)
|
|
112
|
+
Vupper[:-1,1:]*=1/h[:-1]
|
|
113
|
+
V = Vdiag + Vlower + Vupper
|
|
114
|
+
|
|
102
115
|
V[-1][0] = h[-1]
|
|
103
116
|
V[0][-1] = h[-1]
|
|
104
117
|
|
|
@@ -138,5 +151,7 @@
|
|
|
138
151
|
plt.plot(x, y, marker='o', linestyle='None', label='original points')
|
|
139
152
|
plt.legend()
|
|
140
153
|
plt.savefig("graph.png", dpi=200,figsize=(4,3))
|
|
154
|
+
|
|
155
|
+
|
|
141
156
|
```
|
|
142
157
|
|
1
閉曲線verも追加
answer
CHANGED
|
@@ -62,4 +62,81 @@
|
|
|
62
62
|
plt.legend()
|
|
63
63
|
plt.savefig("graph.png", dpi=200,figsize=(4,3))
|
|
64
64
|
```
|
|
65
|
-
|
|
65
|
+
|
|
66
|
+
|
|
67
|
+
周期3次スプライン曲線を使って閉曲線にしました。
|
|
68
|
+
[Periodic cubic smoothing splines as a quadratic
|
|
69
|
+
minimization problem](http://massimozanetti.altervista.org/files/mydocs/periodicCubicSmoothSplines.pdf)
|
|
70
|
+
```python
|
|
71
|
+
import numpy as np
|
|
72
|
+
import matplotlib.pyplot as plt
|
|
73
|
+
|
|
74
|
+
# 周期3次スプライン曲線のパラメータを求める。
|
|
75
|
+
def calculatePeriodicSplineParameter(x):
|
|
76
|
+
# t = 0,1,2,...,x.size-1とする。
|
|
77
|
+
n = x.size-1
|
|
78
|
+
h = np.ones(n)
|
|
79
|
+
a = x
|
|
80
|
+
|
|
81
|
+
# 対角成分
|
|
82
|
+
Adiag = np.eye(n)
|
|
83
|
+
Adiag[1:,1:] *= 2*(h[:-1]+h[1:])
|
|
84
|
+
Adiag[0][0] *= 2*(h[-1]+h[0])
|
|
85
|
+
|
|
86
|
+
# 対角成分の一つ下の斜めの成分
|
|
87
|
+
Alower = np.eye(n, k=-1)
|
|
88
|
+
Alower[1:,:-1]*=h[:-1]
|
|
89
|
+
|
|
90
|
+
# 対角成分の一つ上の斜めの成分
|
|
91
|
+
Aupper = np.eye(n, k=1)
|
|
92
|
+
Aupper[:-1,1:]*=h[:-1]
|
|
93
|
+
A = Adiag + Alower + Aupper
|
|
94
|
+
|
|
95
|
+
A[-1][0] = h[-1]
|
|
96
|
+
A[0][-1] = h[-1]
|
|
97
|
+
|
|
98
|
+
# Aの逆行列
|
|
99
|
+
Ainv = np.linalg.inv(A)
|
|
100
|
+
|
|
101
|
+
V = np.eye(n)*(-2)+np.eye(n, k=1)+np.eye(n, k=-1)
|
|
102
|
+
V[-1][0] = h[-1]
|
|
103
|
+
V[0][-1] = h[-1]
|
|
104
|
+
|
|
105
|
+
r = 3*np.dot(V,a[:-1])
|
|
106
|
+
|
|
107
|
+
c = np.dot(Ainv, r)
|
|
108
|
+
|
|
109
|
+
c = np.insert(c, c.size, c[0])
|
|
110
|
+
d = (c[1:]-c[:-1])/(3*h)
|
|
111
|
+
b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
|
|
112
|
+
a = a[:-1]
|
|
113
|
+
c = c[:-1]
|
|
114
|
+
return [a,b,c,d]
|
|
115
|
+
|
|
116
|
+
# 三次関数の係数から, 元の点から間の点を求める。
|
|
117
|
+
def calculateInterpolationPoints(splineParameter, interpolationPointsNumber=1000):
|
|
118
|
+
a, b, c, d = splineParameter
|
|
119
|
+
assert(a.size==b.size==c.size==d.size)
|
|
120
|
+
index = np.repeat(np.arange(a.size), interpolationPointsNumber, axis=0)
|
|
121
|
+
t = np.tile(np.linspace(0, 1, interpolationPointsNumber), a.size)
|
|
122
|
+
return a[index] + b[index]*t + c[index]*t**2 + d[index]*t**3
|
|
123
|
+
|
|
124
|
+
# 補間したい点の座標
|
|
125
|
+
x = np.random.rand(10)
|
|
126
|
+
y = np.random.rand(x.size)
|
|
127
|
+
x = np.insert(x,x.size,x[0])
|
|
128
|
+
y = np.insert(y,y.size,y[0])
|
|
129
|
+
|
|
130
|
+
parameterX = calculatePeriodicSplineParameter(x)
|
|
131
|
+
parameterY = calculatePeriodicSplineParameter(y)
|
|
132
|
+
|
|
133
|
+
|
|
134
|
+
X = calculateInterpolationPoints(parameterX)
|
|
135
|
+
Y = calculateInterpolationPoints(parameterY)
|
|
136
|
+
|
|
137
|
+
plt.plot(X, Y, label="spline curve")
|
|
138
|
+
plt.plot(x, y, marker='o', linestyle='None', label='original points')
|
|
139
|
+
plt.legend()
|
|
140
|
+
plt.savefig("graph.png", dpi=200,figsize=(4,3))
|
|
141
|
+
```
|
|
142
|
+
|