ベクトルの内積関連の問題を Julia の SymPy で解く

続編書きました
https://blog.goo.ne.jp/r-de-r/e/69e78b604d5a75003c716cd20c47e9de

 

高校数学の「ベクトルの内積」関連の問題をPythonで解く

https://qiita.com/code0327/items/b8614649757110c654b6

「 高校数学の「ベクトルの内積」関連の問題をPythonで解く 」を参考に、値を最後に代入してみた。
https://qiita.com/mrrclb48z/items/2d67795a90778b1ea2d3

Julia の SymPy でやってみる。

問題 1. 次のベクトル a, b の内積と，それらがなす 角θ (0 ≦ θ ≦ π) を求めよ。

(1) a = [2, -3], b = [-4, 6]

まずは，ベクトルの設定からであるが，Python のように sympyMatrix は使わない。公式マニュアルにも，"Avoid sympy.Matrix" と書いてある。

Sym[a, b, c ...] で，Sym 型の配列（今の場合はベクトル）が定義される。

注意すべきは，Sym は関数ではないということ。例えば，Julia では，Int[2, 3] と書くと，2 と 3 を要素とする整数ベクトルが定義できるが，それと同じである。Sym は [a, b, c, ...] の型を示しているのである。

なお，SymPy では（Python の場合も同じであるが），変数の宣言が必要なのは右辺に初めてで来る場合であり，左辺で初めてでてくの変数は宣言する必要はない（しても不都合はないが）。つまり，右辺の変数は左辺の変数の型を持つ変数として宣言されると同時に定義もされるわけである。

  using SymPy
  a = Sym[2, -3];

左辺に初めて出てくる a は，Sym 型の要素を持つベクトル変数として宣言されると同時に定義される。

REPL では自動的に LaTeX 形式で表示されるが，文字型で表示するために "オブジェクト |> string" と言う形式で指示する。

  a |> string

  "Sym[2, -3]"

a は Sym 型のベクトル（1 次元配列）である。

  typeof(a)

  Vector{Sym} (alias for Array{Sym, 1})

  b = Sym[-4, 6]
  b |> string

  "Sym[-4, 6]"

関数などを使用する場合，Python の sympy では，場合によって，どのパッケージ中の関数であるかを明示する必要があるが，Julia の SymPy では，式の評価は対象とするオブジェクトの型によって自動的に適切に行われる。

a, b の内積は a'b，a のノルムは sqrt(a'a) で求めることができる。a'b は a' * b ≡ transpose(a) * b と同じである。

  a'a |> string # a' * a と同じ

  "25"

  sqrt(a'a) |> string # ノルム

  "5"

  a'b |> string # 内積

  "sqrt(x) + 7"

(1) の場合は直接以下によって解 π が求められる。

  acos((a'b)./(sqrt(a'a)*sqrt(b'b))) |> string

  "pi"

(2) a = [1, 1], b = [1-√3, 1+√3]

  using SymPy
  a = Sym[1, 1]
  a |> string

  "Sym[1, 1]"

b の定義で直接 sqrt(3) を使うと，数値解を求めようとするので，それを避けるために変数 x を用いて sqrt(x) とする。x は右辺で初めて出てくるので　@syms で宣言する。

  @syms x
  b = Sym[1-sqrt(x), 1+sqrt(x)]
  b |> string

  "Sym[1 - sqrt(x), sqrt(x) + 1]"

説明のために，段階的に数式を展開していく。

  acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b))) |> string

  "acos(sqrt(2)/sqrt((1 - sqrt(x))*(1 - conjugate(sqrt(x))) + (sqrt(x) + 1)*(conjugate(sqrt(x)) + 1)))"

この段階ではまだ sqrt(x) が実数なのか虚数なのかわかっていないので式が複雑になっている。 そこで，x に 3 を代入する。

  acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b)).(x => 3)) |> string

  "acos(sqrt(2)/sqrt((1 - sqrt(3))^2 + (1 + sqrt(3))^2))"

まだ，簡約化されていないので，更に simplyfi() すると，きれいな結果になる。

  simplify(acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b)).(x => 3))) |> string # "pi/3"

  "pi/3"

もし最初から sqrt(3) でやってしまうと，数値解を求めようとするし，中途で評価が止まってしまう。

  a = Sym[1, 1]
  b = Sym[1-sqrt(3), 1+sqrt(3)]
  acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b))) |> string

  "acos(0.353553390593274*sqrt(2))"

acos の中は

  0.353553390593274*sqrt(2)

  0.5000000000000003

直接 0.5 を与えると，数値解が標示されてしまう。

  acos(0.5)

  1.0471975511965979

この数値が π/3 であることは，知っている人は知っている。

  pi/3

  1.0471975511965976

問題 2. 次の ２ つのベクトルが張る平行四辺形の面積を求めよ。

a = [1, 1], b = [1 - sqrt(3), 1 + sqrt(3)]

平行四辺形の面積は，sqrt(a'a * b'b - (a'b)^2) で求められるが，問題 1. の場合と同じく，いきなり sqrt(3) としてしまうと数値解の方へ流れてしまう。

  using SymPy
  @syms x
  a = Sym[1, 1]
  b = Sym[1-sqrt(x), 1+sqrt(x)]
  b |> string

  "Sym[1 - sqrt(x), sqrt(x) + 1]"

  sqrt(a'a * b'b - (a'b)^2) |> string

  "sqrt(2*(1 - sqrt(x))*(1 - conjugate(sqrt(x))) + 2*(sqrt(x) + 1)*(conjugate(sqrt(x)) + 1) - 4)"

  sqrt(a'a * b'b - (a'b)^2).(x => 3) |> string

  "sqrt(-4 + 2*(1 - sqrt(3))^2 + 2*(1 + sqrt(3))^2)"

  simplify(sqrt(a'a * b'b - (a'b)^2).(x => 3)) |> string

  "2*sqrt(3)"

問題 3. 次のベクトル a に垂直な単位ベクトル ｕ を求めよ。

a = [3, 4]

u = [x, y] とし，
a ⊥ y なので，a'u = 0　 ……　(1)
|u| = 1 なので，sqrt(ｘ^2 + ｙ^2) = 1 　……　(2)
の連立方程式を解く。

  using SymPy
  @syms x, y
  a = Sym[3, 4]
  u = Sym[x, y];

  unorm = sqrt(x^2 + y^2)
  unorm |> string

  "sqrt(x^2 + y^2)"

  answer = solve([a'u, unorm - 1], [x, y])
  answer |> string

  "Tuple{Sym, Sym}[(-4/5, 3/5), (4/5, -3/5)]"

問題 4. 次のベクトル a を法線ベクトルとして，点 (3,2) を通る直線の方程式を求めよ。

問題 3. に引き続き，answer を引用して，

  c1 = answer[1][1]
  c2 = answer[1][2]
  @syms x, y
  solve((x - 3) / c1 - (y - 2) / c2, y) |> string

  "Sym[17/4 - 3*x/4]"

まとめ

Julia の SymPy のほうが簡潔に書けるようだ。