2017年5月20日土曜日

R言語で連立方程式を解く

R言語には連立方程式を解く関数があった

最近R言語にはまってる。
行列を行列として扱える言語のため、逆行列を求める組み込み関数とかあり
逆行列から連立方程式を求めるプログラムが簡単に作れるやん!
って思ってたら、連立方程式を求める組み込み関数が存在してた。
以下の連立方程式で使用方法紹介 \[ \left\{ \begin{array} _2x_1 + 4x_2 = 10 & \\ x_1 + 3x_2 = 6 & \end{array} \right. \tag{1} \]
> # 係数行列を生成
> A <- matrix(c(2, 1, 4, 3),2,2)
> A # 確認
     [,1] [,2]
[1,]    2    4
[2,]    1    3
>
> # 右辺の行列を生成
> b <- matrix(c(10, 6), 2, 1)
> b # 確認
     [,1]
[1,]   10
[2,]    6
>
> #連立方程式を解く組み込み関数
> solve(A,b)
     [,1]
[1,]    3
[2,]    1

超楽で速い!!

ちなみに、この関数的には第二引数(右辺の行列)は配列でもよい 個人的に、行列の連立方程式を意識したかったため、上では第二引数を行列にしてるだけ
> # 配列Ver
>
> # 右辺の行列を配列で生成
> b <- c(10, 6)
> b # 確認
[1] 10  6
>
> #連立方程式を解く組み込み関数
> solve(A,b)
[1] 3 1

使うならこっちかな

拡大係数行列を引数に解を出力する関数

組み込み関数があるのば便利だが、拡大係数行列を引数に解きたい
なので関数にする
仕様としては拡大係数行列を以下のCSV形式で読み込み解を導く感じ
2, 4, 10,
1, 3, 6

まず、この情報から何元の連立方程式か求める必要がある $N$元の場合この配列の要素数は$N(N+1)$なので要素数$L$とすると \[ N(N+1)=L \\ N^2 + N -L = 0 \\ N = \frac{-1+\sqrt{1+4L}}{2} \]
※$L$は正の数となるため
あとは行列制御で何とかなる
#
# 連立方程式を求める関数
#
func <- function(...) {
    data <- c(...)
    # 元の数を求める
    n <- (-1 + sqrt(1 + 4 * length(data))) / 2

    # 拡大係数行列を生成
    M <- t(matrix(data, n+1, n))

    # 係数行列(拡大係数行列から末列を消す)
    A <- M[,-(n+1)]

    # 右辺の行列
    b <- M[,n+1]

    return (solve(A,b))
}

実装も楽すぎて笑うw

試しに13元連立方程式
> func(
+ 24, 67, 77, 7, -4, 43, 72, 22, 61, -5, 52, 24, -19, 2447,
+ 1, 35, 69, 79, -1, 55, 2, 62, 29, 19, 55, -10, 66, 3223,
+ 21, -8, -10, 0, 60, 43, 55, 65, 44, 60, 22, -14, 75, 3483,
+ -6, 59, 28, 73, 72, -6, -13, 11, 69, 6, -16, 76, 2, 2252,
+ 37, -7, 68, 73, 3, 42, 53, -7, 1, 43, 19, 36, 1, 2194,
+ 47, -13, 33, 58, 61, 13, 55, -13, 46, 78, 17, 0, -19, 2150,
+ 23, 48, 33, 64, 6, 6, 76, 6, -19, 56, -13, 59, -17, 1853,
+ 4, 0, 58, 15, 33, 16, 3, 24, 27, 33, -4, -19, -10, 883,
+ 67, 24, 48, 18, 50, 21, 37, 78, 4, -8, 56, 43, 10, 2808,
+ 18, 28, 55, 10, 32, 77, 63, -14, 64, 69, 37, 0, 23, 3202,
+ 75, 77, 69, 21, 27, 36, 54, 78, 45, 74, 12, 4, 73, 4147,
+ 68, -14, 33, 54, 46, -18, 70, 27, 37, -12, 6, -4, 32, 1830,
+ 42, -10, 29, 52, 35, 70, 70, 24, 1, 77, 9, 37, 41, 3449
+ )
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13

できた!!

0 件のコメント:

コメントを投稿