編集履歴

回答編集履歴

整形

2020/07/26 07:05

投稿

スコア698

test CHANGED Viewed

@@ -144,13 +144,15 @@
 # 周期3次スプライン曲線のパラメータを求める。
 def calculatePeriodicSplineParameter(x):
     # t = 0,1,2,...,x.size-1とする。
-    n = x.size-1
+    n = x.size - 1
     h = np.ones(n)
@@ -160,125 +162,103 @@
     # 対角成分
-    Adiag = np.eye(n)
+    Adiag = np.zeros((n, n))
-    Adiag[1:,1:] *= 2*(h[:-1]+h[1:])
+    Adiag[1:, 1:] = np.diag(2 * (h[:-1] + h[1:]))
-    Adiag[0][0] *= 2*(h[-1]+h[0])
+    Adiag[0][0] = 2 * (h[-1] + h[0])
-    # 対角成分の一つ下の斜めの成分
-    Alower = np.eye(n, k=-1)
+    Alower = np.diag(h[:-1], k=-1)
-    Alower[1:,:-1]*=h[:-1]
-    # 対角成分の一つ上の斜めの成分
-    Aupper = np.eye(n, k=1)
+    Aupper = np.diag(h[:-1], k=1)
-    Aupper[:-1,1:]*=h[:-1]
     A = Adiag + Alower + Aupper
     A[-1][0] = h[-1]
     A[0][-1] = h[-1]
-    # Aの逆行列
-    Ainv = np.linalg.inv(A)
     # 対角成分
-    Vdiag = np.eye(n)
+    Vdiag = np.zeros((n, n))
-    Vdiag[1:,1:] *= -(1/h[:-1]+1/h[1:])
+    Vdiag[1:, 1:] = np.diag(-(1 / h[:-1] + 1 / h[1:]))
-    Vdiag[0][0] *= -(1/h[-1]+1/h[0])
+    Vdiag[0][0] = -(1 / h[-1] + 1 / h[0])
-    # 対角成分の一つ下の斜めの成分
-    Vlower = np.eye(n, k=-1)
-    Vlower[1:,:-1]*=1/h[:-1]
+    Vlower = np.diag(1 / h[:-1], k=-1)
-    # 対角成分の一つ上の斜めの成分
-    Vupper = np.eye(n, k=1)
-    Vupper[:-1,1:]*=1/h[:-1]
+    Vupper = np.diag(1 / h[:-1], k=1)
     V = Vdiag + Vlower + Vupper
     V[-1][0] = h[-1]
     V[0][-1] = h[-1]
-    r = 3*np.dot(V,a[:-1])
+    r = 3 * np.dot(V, a[:-1])
-    c = np.dot(Ainv, r)
+    c = np.linalg.solve(A, r)
     c = np.insert(c, c.size, c[0])
-    d = (c[1:]-c[:-1])/(3*h)
+    d = (c[1:] - c[:-1]) / (3 * h)
-    b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
+    b = (a[1:] - a[:-1]) / h - h / 3 * (2 * c[:-1] + c[1:])
     a = a[:-1]
     c = c[:-1]
-    return [a,b,c,d]
+    return [a, b, c, d]
 # 三次関数の係数から, 元の点から間の点を求める。
-def calculateInterpolationPoints(splineParameter, interpolationPointsNumber=1000):
+def calculateInterpolationPoints(splineParameter,
+                                 interpolationPointsNumber=1000):
     a, b, c, d = splineParameter
-    assert(a.size==b.size==c.size==d.size)
+    assert (a.size == b.size == c.size == d.size)
     index = np.repeat(np.arange(a.size), interpolationPointsNumber, axis=0)
     t = np.tile(np.linspace(0, 1, interpolationPointsNumber), a.size)
-    return a[index] + b[index]*t + c[index]*t**2 + d[index]*t**3
+    return a[index] + b[index] * t + c[index] * t**2 + d[index] * t**3
 # 補間したい点の座標
-x = np.random.rand(10)
+x = np.random.rand(100)
 y = np.random.rand(x.size)
-x = np.insert(x,x.size,x[0])
+x = np.insert(x, x.size, x[0])
-y = np.insert(y,y.size,y[0])
+y = np.insert(y, y.size, y[0])
@@ -288,8 +268,6 @@
 X = calculateInterpolationPoints(parameterX)
 Y = calculateInterpolationPoints(parameterY)
@@ -302,11 +280,9 @@
 plt.legend()
-plt.savefig("graph.png", dpi=200,figsize=(4,3))
+plt.savefig("graph.png", dpi=200, figsize=(4, 3))
 ```

コード修正

2020/07/26 07:05

投稿

Penpen7

スコア698

test CHANGED Viewed

@@ -198,7 +198,33 @@
+    # 対角成分
+    Vdiag = np.eye(n)
+    Vdiag[1:,1:] *= -(1/h[:-1]+1/h[1:])
+    Vdiag[0][0] *= -(1/h[-1]+1/h[0])
+    # 対角成分の一つ下の斜めの成分
+    Vlower = np.eye(n, k=-1)
+    Vlower[1:,:-1]*=1/h[:-1]
+    # 対角成分の一つ上の斜めの成分
-    V = np.eye(n)*(-2)+np.eye(n, k=1)+np.eye(n, k=-1)
+    Vupper = np.eye(n, k=1)
+    Vupper[:-1,1:]*=1/h[:-1]
+    V = Vdiag + Vlower + Vupper
     V[-1][0] = h[-1]
@@ -278,6 +304,10 @@
 plt.savefig("graph.png", dpi=200,figsize=(4,3))
 ```
 ![イメージ説明](b25297d263933355338c31ffe468e743.png)

閉曲線verも追加

2020/07/26 05:56

投稿

Penpen7

スコア698

test CHANGED Viewed

@@ -127,3 +127,157 @@
 ```
 ![イメージ説明](97c19caa7151c4bbf2485e2189f77066.png)
+周期3次スプライン曲線を使って閉曲線にしました。
+[Periodic cubic smoothing splines as a quadratic
+minimization problem](http://massimozanetti.altervista.org/files/mydocs/periodicCubicSmoothSplines.pdf)
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+# 周期3次スプライン曲線のパラメータを求める。
+def calculatePeriodicSplineParameter(x):
+    # t = 0,1,2,...,x.size-1とする。
+    n = x.size-1
+    h = np.ones(n)
+    a = x
+    # 対角成分
+    Adiag = np.eye(n)
+    Adiag[1:,1:] *= 2*(h[:-1]+h[1:])
+    Adiag[0][0] *= 2*(h[-1]+h[0])
+    # 対角成分の一つ下の斜めの成分
+    Alower = np.eye(n, k=-1)
+    Alower[1:,:-1]*=h[:-1]
+    # 対角成分の一つ上の斜めの成分
+    Aupper = np.eye(n, k=1)
+    Aupper[:-1,1:]*=h[:-1]
+    A = Adiag + Alower + Aupper
+    A[-1][0] = h[-1]
+    A[0][-1] = h[-1]
+    # Aの逆行列
+    Ainv = np.linalg.inv(A)
+    V = np.eye(n)*(-2)+np.eye(n, k=1)+np.eye(n, k=-1)
+    V[-1][0] = h[-1]
+    V[0][-1] = h[-1]
+    r = 3*np.dot(V,a[:-1])
+    c = np.dot(Ainv, r)
+    c = np.insert(c, c.size, c[0])
+    d = (c[1:]-c[:-1])/(3*h)
+    b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
+    a = a[:-1]
+    c = c[:-1]
+    return [a,b,c,d]
+# 三次関数の係数から, 元の点から間の点を求める。
+def calculateInterpolationPoints(splineParameter, interpolationPointsNumber=1000):
+    a, b, c, d = splineParameter
+    assert(a.size==b.size==c.size==d.size)
+    index = np.repeat(np.arange(a.size), interpolationPointsNumber, axis=0)
+    t = np.tile(np.linspace(0, 1, interpolationPointsNumber), a.size)
+    return a[index] + b[index]*t + c[index]*t**2 + d[index]*t**3
+# 補間したい点の座標
+x = np.random.rand(10)
+y = np.random.rand(x.size)
+x = np.insert(x,x.size,x[0])
+y = np.insert(y,y.size,y[0])
+parameterX = calculatePeriodicSplineParameter(x)
+parameterY = calculatePeriodicSplineParameter(y)
+X = calculateInterpolationPoints(parameterX)
+Y = calculateInterpolationPoints(parameterY)
+plt.plot(X, Y, label="spline curve")
+plt.plot(x, y, marker='o', linestyle='None', label='original points')
+plt.legend()
+plt.savefig("graph.png", dpi=200,figsize=(4,3))
+```
+![イメージ説明](b25297d263933355338c31ffe468e743.png)