質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.90%

Pythonで行列の計算をしたい。

解決済

回答 4

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 286

Pythonで行列の計算をし,パラメータを算出したいです.xとyの値は例として挙げています.予定ではクリックで点を生成し,その点をふやしていくことを考えています.つまり点の個数が決まっていなくても,図のような行列の計算はできますでしょうか.最終的な結果としてcというパラメータの値を取り出したいです.
よろしくお願い致します.

(追記)
手計算によって求めた解を載せます。また,これらの行列は離散点を3次スプライン補間をする際の3次多項式のパラメータの一部(c)を求める式になります。
よろしくお願い致します。

x y
0 0
1 3
4 3
5 1
8 2

イメージ説明

手計算によって求めた解

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • can110

    2020/07/23 18:53 編集

    まずは点の数がたとえば5個と決まっている場合はコードは書けますでしょうか?
    書けない場合でも、A,bおよびcを手計算で求めた結果を提示ください。
    さらに、この問題、数式に名前(~の式、~の解法など)がある場合はそれを提示ください。

    キャンセル

  • can110

    2020/07/24 02:05

    Aの最も右下の値は0ではなく1のはずではないでしょうか?
    またb=[ 0 -9 -6 7 0]ではないでしょうか?
    今一度、数式または手計算の結果が正しいか確認ください。

    キャンセル

  • takes.it.easy

    2020/07/24 12:53

    どちらもおっしゃる通りです。再度計算すると
    b=[0, -9, -6, 7, 0]
    c=[0, -0.930556, -0.518519, 0.939815, 0]
    となりました。

    キャンセル

回答 4

+2

x,yにリストとして座標を追加し、numpyの配列に変換すれば以下のソースコードで計算できると思います。
Ax=rとして解いた。a,b,c,dは三次関数の係数。

import numpy as np

x = np.array([0, 1, 4, 5, 8])
y = np.array([0, 3, 3, 1, 2])
n = len(x)-1
h = x[1:]-x[:-1]
a = y

# 対角成分
Adiag = np.eye(n+1)
Adiag[1:-1, 1:-1] *= 2*(h[:-1]+h[1:])

# 対角成分の一つ下の斜めの成分
Alower = np.eye(n+1, k=-1)
Alower[1:, :-1] *= h
Alower[-1, -2] = 0

# 対角成分の一つ上の斜めの成分
Aupper = np.eye(n+1, k=1)
Aupper[:-1, 1:] *= h
Aupper[0, 1] = 0

A = Adiag + Alower + Aupper

# Aの逆行列
Ainv = np.linalg.inv(A)

r = np.zeros(n+1)
r[1:-1] = (3/h)[1:]*(a[2:]-a[1:-1])-(3/h)[:-1]*(a[1:-1]-a[:-2])

c = np.dot(Ainv, r)

d = (c[1:]-c[:-1])/(3*h)
b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])
print("x", x, "y", y)
print("h", h, "a", a)
print("A")
print(A)
print("r", r)

print("a", a)
print("b", b)
print("c", c)
print("d", d)
x [0 1 4 5 8] y [0 3 3 1 2]
h [1 3 1 3] a [0 3 3 1 2]
A
[[1. 0. 0. 0. 0.]
 [1. 8. 3. 0. 0.]
 [0. 3. 8. 1. 0.]
 [0. 0. 1. 8. 3.]
 [0. 0. 0. 0. 1.]]
r [ 0. -9. -6.  7.  0.]
a [0 3 3 1 2]
b [ 3.31018519  2.37962963 -1.96759259 -1.5462963 ]
c [ 0.         -0.93055556 -0.51851852  0.93981481  0.        ]
d [-0.31018519  0.04578189  0.48611111 -0.10442387]


イメージ説明
イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/24 17:29

    ありがとうございます。
    スプライン曲線を描く際,このコードではノンパラメトリックとして計算されていると思いますが,xは短調増加でなくてはならない制限を解決したく、パラメトリックの式として描画したいと考えています。パラメトリックな3次多項式とするにはどのようにすればよういでしょうか。変数をtとおいたときのx,yの式やパラメータの算出などご教授いただけると幸いです。

    キャンセル

  • 2020/07/24 18:04 編集

    パラメトリックの式ではxは任意の順序でもスプライン曲線を描けるという理由が申し訳ありませんが、わかりませんでした...
    それよりもxが常に単調増加するように座標を挿入するところを選んで挿入するか、スプライン曲線を描く前にxとyの対応を崩さないようにxをソートしてから各パラメータを求めるようにするのが自然なように感じます。例えば, xとyの対応を崩さないようにソートする場合は
    y = y[np.argsort(x)]
    x = np.sort(x)
    でxを単調増加数列に変えることができます。
    なお、スプライン曲線は三次関数を滑らかに繋いでいるにすぎないため、三次関数をパラメータ表示できるかということと同値ですが、それは可能だと思います。(が、やって効果はあるかは不明です)

    キャンセル

  • 2020/07/25 23:04

    媒介変数表示しても一切やることは変わらないので(文字が入れ替わるだけ)、回答を参考に改変してください。

    キャンセル

checkベストアンサー

0

パラメータ表示でやってみた。

import numpy as np
import matplotlib.pyplot as plt

# スプライン曲線のパラメータを求める。
def calculateSplineParameter(x):
    # t = 0,1,2,...,x.size-1とする。
    n = x.size-1
    h = np.ones(n)
    a = x

    # 対角成分
    Adiag = np.eye(n+1)
    Adiag[1:-1, 1:-1] *= 2*(h[:-1]+h[1:])

    # 対角成分の一つ下の斜めの成分
    Alower = np.eye(n+1, k=-1)
    Alower[1:, :-1] *= h
    Alower[-1, -2] = 0

    # 対角成分の一つ上の斜めの成分
    Aupper = np.eye(n+1, k=1)
    Aupper[:-1, 1:] *= h
    Aupper[0, 1] = 0

    A = Adiag + Alower + Aupper

    # Aの逆行列
    Ainv = np.linalg.inv(A)

    r = np.zeros(n+1)
    r[1:-1] = (3/h)[1:]*(a[2:]-a[1:-1])-(3/h)[:-1]*(a[1:-1]-a[:-2])

    c = np.dot(Ainv, r)
    d = (c[1:]-c[:-1])/(3*h)
    b = (a[1:]-a[:-1])/h-h/3*(2*c[:-1]+c[1:])

    return [a,b,c,d]

# 三次関数の係数から, 元の点から間の点を求める。
def calculateInterpolationPoints(splineParameter, orignalPointsNumber, interpolationPointsNumber=1000):
    a, b, c, d = splineParameter

    index = np.repeat(np.arange(orignalPointsNumber-1), interpolationPointsNumber, axis=0)
    t = np.tile(np.linspace(0, 1, interpolationPointsNumber), orignalPointsNumber-1)

    return a[index] + b[index]*t + c[index]*t**2 + d[index]*t**3

# 補間したい点の座標
x = np.random.rand(20)
y = np.random.rand(x.size)

parameterX = calculateSplineParameter(x)
parameterY = calculateSplineParameter(y)

X = calculateInterpolationPoints(parameterX, x.size)
Y = calculateInterpolationPoints(parameterY, y.size)

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))


イメージ説明

周期3次スプライン曲線を使って閉曲線にしました。
Periodic cubic smoothing splines as a quadratic minimization problem

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.zeros((n, n))
    Adiag[1:, 1:] = np.diag(2 * (h[:-1] + h[1:]))
    Adiag[0][0] = 2 * (h[-1] + h[0])
    Alower = np.diag(h[:-1], k=-1)
    Aupper = np.diag(h[:-1], k=1)

    A = Adiag + Alower + Aupper
    A[-1][0] = h[-1]
    A[0][-1] = h[-1]

    # 対角成分
    Vdiag = np.zeros((n, n))
    Vdiag[1:, 1:] = np.diag(-(1 / h[:-1] + 1 / h[1:]))
    Vdiag[0][0] = -(1 / h[-1] + 1 / h[0])
    Vlower = np.diag(1 / h[:-1], k=-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])

    c = np.linalg.solve(A, 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(100)
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))


イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/27 12:16

    ありがとうございます。
    勉強させていただきます。

    キャンセル

0

できるかとという問いでしたら、できます。
numpy.insert()とかで行と列を拡張していけばいいんではないでしょうか。

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

予定ではクリックで点を生成し,その点をふやしていくことを考えています.

上記を実現するため、クラスを使って以前の状態を保持しておく例です。
なお、点の追加において、x値`は単純増加で与える(制御点が左側にぐるっと戻らない)ことが前提です。

import numpy as np
import matplotlib.pyplot as plt

# 3次スプライン補間法
# http://www.civil.kumamoto-u.ac.jp/matsu/spline.pdf

class Spline:
    # xys : 最低3点ある配列
    def __init__(self, xys):
        self.xys = np.array(xys)
        N = self.xys.shape[0]
        assert N >= 3

        self.h = np.array([self.xys[i,0]-self.xys[i-1,0] for i in range(1,N)])
        self.a = self.xys[:,1]

        self.A = np.eye(N)
        for i in range(2,N):
            self.A[i-1,i-2:i+1] = self.get_A_line(i)

        self.r = np.zeros(N)
        for i in range(2,N):
            self.r[i-1] = self.get_r_val(i)

    # 指定行のAの値を計算して返す
    def get_A_line(self, i):
        return [self.h[i-2], 2*(self.h[i-2]+self.h[i-1]), self.h[i-1]]

    # 指定位置のrの値を計算して返す
    def get_r_val(self, i):
        return 3*(self.a[i]-self.a[i-1])/self.h[i-1] - 3*(self.a[i-1]-self.a[i-2])/self.h[i-2]

    # 新たな点を追加
    def add(self,xy):
        self.xys = np.vstack([self.xys, np.array(xy)])
        N = self.xys.shape[0]
        self.h = np.append(self.h, self.xys[N-1,0]-self.xys[N-2,0])
        self.a = self.xys[:,1]

        Anew = np.zeros((N,N))
        Anew[0:N-1,0:N-1] = self.A
        self.A = Anew
        self.A[N-2,N-3:N] = self.get_A_line(N-1)
        self.A[N-1,N-1] = 1

        self.r = np.append(self.r, 0)
        self.r[N-2] = self.get_r_val(N-1)

    # パラメータを算出
    def calc_param(self):
        Ainv = np.linalg.inv(self.A)
        self.c = np.dot(Ainv, self.r)
        N = self.xys.shape[0]
        self.b = np.array([ (self.a[i]-self.a[i-1])/self.h[i-1] - self.h[i-1]/3*(2*self.c[i-1]+self.c[i]) for i in range(1,N)])
        self.d = np.array([ (self.c[i]-self.c[i-1])/(3*self.h[i-1]) for i in range(1,N)])

    # 補間されたy値を返す
    def get_y(self, x):
        xs = self.xys[:,0]
        if x < xs[0] or x > xs[-1]:
            return 0
        elif x == self.xys[0,1]:
            return xs[0]
        elif x == xs[-1]:
            return self.xys[-1,1]

        # 区間を決定
        i = 0
        while True:
            if xs[i] > x:
                break
            i += 1
        i -= 1

        dx = (x-xs[i])
        return self.a[i] + self.b[i]*dx + self.c[i]*dx*dx + self.d[i]*dx*dx*dx


sp = Spline([[0,0],[1,3],[4,3]]) # 初期の3点は最初に与える
sp.add([5,1])
sp.add([8,2])
sp.calc_param()

print(sp.r)# [ 0. -9. -6.  7.  0.]
print(sp.a)# [0 3 3 1 2]
print(sp.b)# [ 3.31018519  2.37962963 -1.96759259 -1.5462963 ]
print(sp.c)# [ 0.         -0.93055556 -0.51851852  0.93981481  0.        ]
print(sp.d)# [-0.31018519  0.04578189  0.48611111 -0.10442387]

# 描画
xs = np.linspace( np.min(sp.xys[:,0]), np.max(sp.xys[:,0]), 100)
ys = np.array([sp.get_y(x) for x in xs])
plt.scatter(sp.xys[:,0], sp.xys[:,1])
plt.plot(xs, ys)
plt.show()


イメージ説明

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/24 16:14

    丁寧な回答ありがとうございます。
    こちらのコードではxは単純増加ということですが、自由な曲線を描くために、媒介変数t(0≦t≦1)を用いて、補間したいと考えています。媒介変数でスプライン補間をする場合、それぞれのパラメータはどのように決めるのでしょうか。

    キャンセル

  • 2020/07/24 19:46

    > こちらのコードではxは単純増加ということですが
    もともとの数式、考えが単純増加を前提としていると思うのですが、違いますでしょうか?

    > 自由な曲線を描くために~
    そのようにできる数式、考えはどのようなものでしょうか?

    キャンセル

  • 2020/07/25 13:58

    もともとの式はxが決まるとyが決まってしまうため、コースを一周するような曲線を描くことができません。そこで変数tを用いて媒介変数表示してみたらどうかと考えています。
    t = [0,1]
    h = t_i+1 - t_i

    x(t) = a_x + b_x * t + c_x * t^2 + d_x * t^3
    a_x,i = x_i
    y(t) = a_y + b_y * t + c_y * t^2 + d_y * t^3
    a_y,i = y_i

    よろしくお願いします。
    参考サイト:https://tajimarobotics.com/cubic-spline-interpolation-parameter/

    キャンセル

  • 2020/07/25 14:12

    なるほど。
    まずは固定の点の数で動くコードを書いてみてはいかがでしょうか?

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.90%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る