裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

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

2021年11月21日 | ブログラミング

ベクトルの内積関連の問題を Julia の SymPy で解く
https://blog.goo.ne.jp/r-de-r/e/cd7d0c0701c2541aa17d5943fe6e482b

の問題 2 の (2) で,「a = [1, 1], b = [1-√3, 1+√3]b の定義で直接 sqrt(3) を使うと,数値解を求めようとするので,それを避けるために変数 x を用いて sqrt(x) とする。」などと書いてしまったが,誤ったことを書いてしまった。

ちゃんとやればちゃんとできる。SymPy 公式サイトの文書をよく読んだら(というか,最初の方に)解決法は書いてあった。

sqrt(3) では,Base の sqrt が呼ばれてしまうので,数値解になってゆく。

ついでに,関数を内包していく書き方ではなく,パイプ |> を使ってみる。

  using SymPy
  u = sqrt(3)
  "u = $u, typeof(u) = $(typeof(u))" |> println
  u = 1.7320508075688772, typeof(u) = Float64

ここは,sympy.sqrt(3) として,Sym オブジェクトの √3 にしなければならない。

  v = sympy.sqrt(3)
  "v = $v, typeof(v) = $(typeof(v))" |> println
  v = sqrt(3), typeof(v) = Sym

ということで,正解は以下のとおりである。

  a = Sym[1, 1]
  b = Sym[1-sympy.sqrt(3), 1+sympy.sqrt(3)]
  "a = $a\nb = $b" |> println
  a = Sym[1, 1]
  b = Sym[1 - sqrt(3), 1 + sqrt(3)]
  acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b))) |> println
  acos(sqrt(2)/sqrt((1 - sqrt(3))^2 + (1 + sqrt(3))^2))

もう少し簡約化するために simplify を使う。

  acos((a'*b)./(sqrt(a'*a)*sqrt(b'*b))) |> simplify |> println
  pi/3

閉じカッコに気をつけなくてもよいので,統一感も必用だが,使えるときにはパイプもよいかな。

なお,b の定義のとき,

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

でもよい。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

微分関連の問題を Julia の SymPy で解く

2021年11月21日 | ブログラミング

高校数学の「微分」関連の問題をPythonで解く

https://qiita.com/code0327/items/97fcacc42ab44402eb35

Julia の SymPy と Plotでやってみる。

  using SymPy
  @syms x
  diff((x^2 - 1) * (2x^2 + x - 3)) |> string
  "2*x*(2*x^2 + x - 3) + (4*x + 1)*(x^2 - 1)"

結果を因数分解した場合

  factor(diff((x^2 - 1) * (2x^2 + x - 3))) |> string
  "(x - 1)*(8*x^2 + 11*x + 1)"

結果を展開した場合

  expand(diff((x^2 - 1) * (2x^2 + x - 3))) |> string
  "8*x^3 + 3*x^2 - 10*x - 1"

  using SymPy
  @syms x
  diff((x + 1)^2 * (2x - 3)^4) |>string
  "8*(x + 1)^2*(2*x - 3)^3 + (2*x - 3)^4*(2*x + 2)"

結果を因数分解した場合

  factor(diff((x + 1)^2 * (2x - 3)^4)) |>string
  "2*(x + 1)*(2*x - 3)^3*(6*x + 1)"

結果を式を展開した場合

  expand(diff((x + 1)^2 * (2x - 3)^4)) |>string
  "96*x^5 - 320*x^4 + 160*x^3 + 360*x^2 - 270*x - 54"
  
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

極限値関連の問題を Julia の SymPy で解く

2021年11月21日 | ブログラミング

高校数学の「極限値」関連の問題をPythonで解く

https://qiita.com/code0327/items/1fd0df574c8c8357af5a

Julia の SymPy と Plotでやってみる。

  using SymPy
  @syms x
  limit((5x^2 - x) / (4 -2x^2), x, -oo) |> string
  "-5/2"

  using SymPy
  @syms x
  limit(sqrt(x^2 + 1) - x, x, oo) |>string
  "0"
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

2次関数関連の問題を Julia の SymPy で解く

2021年11月21日 | ブログラミング

高校数学の「2次関数」関連の問題をPythonで解く

https://qiita.com/code0327/items/66efb7ac56cdea5b0b0d

Julia の SymPy と Plotでやってみる。

  using SymPy
  @syms x
  y = -2x^2 + 4x + 1
  y |> string
  "-2*x^2 + 4*x + 1"

y の導関数 y' を求め,y1 と置く。

  y1 = diff(y)
  y1 |> string
  "4 - 4*x"

y1 = 0 を解き,切線の傾きが 0 になるときの x 座標を求める。

  x0 = solve(y1)[1];

そのときの y の値を求める。

  y0 = y(x => x0);
  println("頂点の座標は ($x0, $y0)")
  頂点の座標は (1, 3)

関数の定義

  f(x) = -2x^2 + 4x + 1
  f (generic function with 1 method)

必要最小限の描画

  using Plots
  plot(f, xlims=(-1, 3), label="")
  scatter!([x0], [y0], label="")

  using Plots
  f(x) = x^2 +x - 6
  plot(f, label="")

x 軸との交点を求めて...などしなくて,x^2 + x - 6 < 0 を解くだけでよい。

  using SymPy
  @syms x
  answer = solve(Lt(x^2 + x - 6, 0))
  answer |> string
  "(-3 < x) & (x < 2)"

  f(x) = x^2 + x + 1
  plot(f, label="")

問 2 に同じだが,解がないときは False が返される。

  solve(Lt(x^2 + x + 1, 0)) |> string
  "False"
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

高次方程式・恒等式関連の問題を Julia の SymPy で解く

2021年11月21日 | ブログラミング

1. 高校数学の「高次方程式・恒等式」関連の問題をPythonで解く

https://qiita.com/code0327/items/66efb7ac56cdea5b0b0d

Julia の SymPy でやってみる。


  
  using SymPy
  @syms x
  f = x^3 - 1;
  # (1)
  f.(x => 2) |> string
  "7"
  # (2)
  f.(x => 1) |> string
  "0"
  # (3)
  factor(f) |> string
  "(x - 1)*(x^2 + x + 1)"
  # (4)
  solve(f) |> string
  "Sym[1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]"
  # (1)
  solve(x^4 - 10x^3 +35x^2 - 50x + 24) |> string
  "Sym[1, 2, 3, 4]"
  # なお,
  factor(x^4 - 10x^3 +35x^2 - 50x + 24) |> string
  "(x - 4)*(x - 3)*(x - 2)*(x - 1)"
  # (2)
  solve(x^4 - 10x^2 +9) |> string
  "Sym[-3, -1, 1, 3]"
  # なお,
  factor(x^4 - 10x^2 +9) |> string
  "(x - 3)*(x - 1)*(x + 1)*(x + 3)"
  # (1)
  @syms a, b, c, x
  solve(Eq(x^2, a*(x - 2)^2 + b*(x - 2) + c), [a, b, c]) |> string
  "Dict{Any, Any}(c => 4, b => 4, a => 1)"
  # (2)
  @syms a, b, x
  solve(Eq((3x- 4) / ((x + 2) * (x - 3)), a / (x + 2) + b / (x - 3)), [a, b])
  Dict{Any, Any} with 2 entries:
    b => 1
    a => 2
  # (1)
  solve(Eq(3 / (x + 2) + 3 / (x - 2), 2)) |> string
  "Sym[-1, 4]"
  # (2)
  solve(Eq(3 / (x + 2) + 3 / (x - 2), 0)) |> string
  "Sym[0]"
  # (3)
  solve(Eq(sqrt(x + 3), x + 1)) |> string
  "Sym[1]"
  # (4)
  solve(Eq(-sqrt(x + 3), x + 1)) |> string
  "Sym[-2]"
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

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

2021年11月21日 | ブログラミング

続編書きました
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'ba' * 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. 次の 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 に垂直な単位ベクトル u を求めよ。

a = [3, 4]

u = [x, y] とし,
a ⊥ y なので,a'u = 0  …… (1)
|u| = 1 なので,sqrt(x^2 + y^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 のほうが簡潔に書けるようだ。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村