🦋

Juliaで連立方程式を解く

2023/07/17に公開

Julia

Juliaで連立方程式を解く

はじめに

2023年7月10日〜7月12日に九州大学ので開催された「数学と物理におけるJuliaの活用」に参加しました。7月11日の1日だけの参加（東京→博多→東京の日帰り！）だったのですが，日本でJuliaを薦めている多くの方とお会いできてとても楽しかったです。

Webサイトもありますので，これからJuliaを始める方はとても参考になると思います。

東京に帰ってきて，「何かJuliaでやりたいな。」と思っていたところ，うちの長男（大学1年生）が線形代数の問題に取り組んでいたので，その問題をjuliaで解いてみることにしました。
同じ問題でも良かったのですが，数字は少し変えました。

連立方程式を解く

次の2つの連立方程式を解くことを考えます。

\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-4\\ \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}1\\1\\1\\1\\ \end{pmatrix}

\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-3\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}

1つ目の連立方程式

\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-4\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}1\\1\\1\\1\\ \end{pmatrix}

まず，

A_1=\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-4\\ \end{pmatrix}\qquad B_1=\begin{pmatrix}1\\1\\1\\1\\ \end{pmatrix}

を定義します。

X_1=\begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}

はJuliaではA₁\B₁で求めることができます。

4×4 Matrix{Int64}:
1 2 1 -1
1 2 -1 1
2 1 1 2
1 3 2 -4

4-element Vector{Int64}:
1
1
1
1

A₁\B₁

4-element Vector{Float64}:
-0.3333333333333328
0.6666666666666664
0.3333333333333331
0.3333333333333331

小数表示でもいいのですが，有理数の分数表示をしてみます。

A₁ = [
   1 2 1 -1
   1 2 -1 1
   2 1 1 2
   1 3 2 -4//1
]

4×4 Matrix{Rational{Int64}}:
1//1 2//1 1//1 -1//1
1//1 2//1 -1//1 1//1
2//1 1//1 1//1 2//1
1//1 3//1 2//1 -4//1

B₁ = [
   1
   1
   1
   1//1
]

4-element Vector{Rational{Int64}}:
1//1
1//1
1//1
1//1

A₁\B₁

4-element Vector{Rational{Int64}}:
-1//3
2//3
1//3
1//3

$X_1=\begin{pmatrix}-1/3\\2/3\\1/3\\1/3\\ \end{pmatrix}$ であることがわかりました。

2つ目の連立方程式

\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-3\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}

同じように，

A_2=\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-3\\ \end{pmatrix}\quad B_2=\begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}

を定義します。

X_2=\begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}

はA₂\B₂で求めてみます。

4×4 Matrix{Int64}:
1 2 1 -1
1 2 -1 1
2 1 1 2
1 3 2 -3

B₂= [
   1
   3
   -3
   2
]

4-element Vector{Int64}:
1
3
-3
2

A₂\B₂

4-element Vector{Float64}:
1.6000000000000005
0.19999999999999946
-2.8
-1.8000000000000003

こちらも有理数の分数表示をしてみましょう。

A₂= [
   1 2 1 -1
   1 2 -1 1
   2 1 1 2
   1 3 2 -3//1
]

4×4 Matrix{Rational{Int64}}:
1//1 2//1 1//1 -1//1
1//1 2//1 -1//1 1//1
2//1 1//1 1//1 2//1
1//1 3//1 2//1 -3//1

B₂= [
   1
   3
   -3
   2//1
]

4-element Vector{Rational{Int64}}:
1//1
3//1
-3//1
2//1

A₂\B₂

SingularException(4)

Stacktrace:
[1] checknonsingular
@ ~/.julia/juliaup/julia-1.9.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:19 [inlined]
[2] generic_lufact!(A::Matrix{Rational{Int64}}, pivot::RowMaximum; check::Bool)

......

あれ，エラーだ。どうしてだろう？（わざとらしい。。。）

ということで，ランクを調べてみます。行列のランクを調べるにはLinearAlgebra.jlのパッケージを利用します。

行列のランクを調べ，逆行列を確認

using LinearAlgebra
rank(A₂)

$rank(A_2)=3$ だ。逆行列持たない。どうして，実数の場合のA₂\B₂が求まったのか不思議ですが，恐らく，解が不定で，解の1つが表示されたと思われます。

いずれにして，ちゃんと求めなくては！

rank(A₁)

$rank(A_1)=4$ です。こちらは，ちゃんと逆行列がありますね。一応，有理数バージョンで確認しておきます。

逆行列はinv(A)で求めます。

inv(A₁)

4×4 Matrix{Rational{Int64}}:
-9//2 5//6 4//3 2//1
5//2 -1//6 -2//3 -1//1
5//2 -5//6 -1//3 -1//1
2//1 -1//3 -1//3 -1//1

inv(A₂)

SingularException(4)
.......

$A_1$ の逆行列は求まり， $A_2$ の逆行列は求まらないですね。

行簡約階段形（Reduced row echelon form）

行簡約階段形を行います。詳細は

をみてください。RowEchelon.jlパッケージのrref(A)コマンドを用います。

https://github.com/blegat/RowEchelon.jl

using RowEchelon
rref(A₂)

4×4 Matrix{Rational{Int64}}:
1//1 0//1 0//1 2//1
0//1 1//1 0//1 -1//1
0//1 0//1 1//1 -1//1
0//1 0//1 0//1 0//1

ちょっと目がチカチカするのでLatexify.jlを＼パッケージを用いて，数式っぽくします。

using Latexify
rref(A₂)|>latexify

第4行に0が並び，ランクが3であることがわかります。

まず， $A_2$ の右側に単位行列を並べます。

A₂E = hcat(A₂,Matrix(I, 4 ,4))

4×8 Matrix{Rational{Int64}}:
1//1 2//1 1//1 -1//1 1//1 0//1 0//1 0//1
1//1 2//1 -1//1 1//1 0//1 1//1 0//1 0//1
2//1 1//1 1//1 2//1 0//1 0//1 1//1 0//1
1//1 3//1 2//1 -3//1 0//1 0//1 0//1 1//1

行基本変形を用いて，行簡約階段形にします。

rref(A₂E)

4×8 Matrix{Rational{Int64}}:
1//1 0//1 0//1 2//1 0//1 1//12 7//12 -1//4
0//1 1//1 0//1 -1//1 0//1 1//4 -1//4 1//4
0//1 0//1 1//1 -1//1 0//1 -5//12 1//12 1//4
0//1 0//1 0//1 0//1 1//1 -1//6 -1//6 -1//2

\begin{pmatrix}1& 2& 1& -1 \\ 1 &2 &-1 &1\\ 2 &1 &1 &2\\ 1 &3 &2 &-3\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}

を解くには

\begin{pmatrix}1& 0& 0& 2 \\ 0 &1 &0 &-1\\ 0 &0 &1 &-1\\ 0 &0 &0 &0\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}0& 1/12& 7/12& -1/4 \\ 0 &1/4 &-1/4 &1/4\\ 0 &-5/12 &1/12 &1/4\\ 1 &-1/6 &-1/6 &-1/2\\ \end{pmatrix} \begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}

を解けばOKです。

X_2=\begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}

とおきます。

Symbolics.jlパッケージを用いて，文字を扱えるようにします。

using Symbolics
@variables x₁ x₂ x₃ x₄
X₂=[x₁ 
    x₂ 
    x₃ 
    x₄]

rref(A₂)*X₂

\begin{pmatrix}1& 0& 0& 2 \\ 0 &1 &0 &-1\\ 0 &0 &1 &-1\\ 0 &0 &0 &0\\ \end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}= \begin{pmatrix}x_1+2x_4\\x_2-x_4\\x_3-x_4\\0\\ \end{pmatrix}

となります。

rref(A₂E)[1:4,5:8]*B₂

4-element Vector{Rational{Int64}}:
-2//1
2//1
-1//1
0//1

\begin{pmatrix}0& 1/12& 7/12& -1/4 \\ 0 &1/4 &-1/4 &1/4\\ 0 &-5/12 &1/12 &1/4\\ 1 &-1/6 &-1/6 &-1/2\\ \end{pmatrix} \begin{pmatrix}1\\3\\-3\\2\\ \end{pmatrix}= \begin{pmatrix}-2\\2\\-1\\0\\ \end{pmatrix}

となります。

第4行が共に0であるので， $x_4$ は任意の実数として， $x_1$ ， $x_2$ ， $x_3$ を表すことにすると，

eqs = [(rref(A₂)*X₂)[i]~(rref(A₂E)[1:4,5:8]*B₂)[i] for i =1:3]
var = [x₁, x₂ ,x₃]
Symbolics.solve_for(eqs,var)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
(-2//1) - (2//1)*x₄
(2//1) + x₄
x₄ - (1//1)

X_2=\begin{pmatrix}x_1\\x_2\\x_3\\x_4\\ \end{pmatrix}=\begin{pmatrix}-2-2x_4\\2-x_4\\x_4-1\\x_4\\ \end{pmatrix}

と求まりました。

$x_4=-1.8$ とすると， $x_1=1.6$ ， $x_2=0.2$ ， $x_3=-2.8$ となり，最初に実数でを使って求めた値と一致します。

まとめ

連立方程式を解くときは，まずランクを調べよう！

Juliaで連立方程式を解く

はじめに

連立方程式を解く

1つ目の連立方程式

2つ目の連立方程式

行列のランクを調べ，逆行列を確認

行簡約階段形（Reduced row echelon form）

まとめ

Discussion