2 直線の交点 SymPy/Julia

あいも変わらず，

2直線の交点座標を求めるプログラム
https://blog.goo.ne.jp/r-de-r/e/e81316c26c521b31e9d47526f9bd5861

がよく参照されているが，わけわからん。

そこで，Julia の SymPy で書いたものも掲載しておく。
これが「わけわからん」という人もいるだろうが。

using SymPy

@syms a b c d x

function func(a, b, c, d)
    if a != c
        eq1 = a*x + b
        eq2 = c*x + d
        x2 = solve(Eq(eq1, eq2), x)[1]
        y2 = simplify(eq1(x => x2))
        println("x = $x2,  y = $y2")
    elseif b == d
        println("2直線はまったく同じなので，交点は無数にあります")
    else
        println("2直線は平行なので，交点はありません")
    end
end

記号変数で呼び出すと，数式で答えが返ってくる。

func(a, b, c, d)
x = (-b + d)/(a - c),  y = (a*d - b*c)/(a - c)

数値で呼び出すと，然るべき答えが返ってくる。

func(2, 3, 1, -4)
x = -7,  y = -11

func(1, 2, 1, 2)
2直線はまったく同じなので，交点は無数にあります

func(1, 2, 1, 5)
2直線は平行なので，交点はありません

SymPy/Julia(2)

SymPy/Julia

https://www.math.csi.cuny.edu/~verzani/tmp/julia/symbolic.html

5.　数式の図示

数式は，Julia へ戻ることなく，直接プロットすることができる。

f(x) = sin(x^4) - cos(x^4) を [0,1] の範囲で描くには以下のようにする。

  using Plots
  gr(label="")
  f(x) = sin(x^4) - cos(x^4)
  plot(f(x), 0, 1)

plot の第一引数は関数でなくて，数式でよい。

  g = sin(x^4) - cos(x^4)
  plot(g, 0, 1)

f(x) も g も，どちらも Sym オブジェクトである。

  println("f(x): $(f(x)),  typeof(f(x)): $(typeof(f(x)))")
  println("g:    $g,  typeof(g):    $(typeof(g))")

  f(x): sin(x^4) - cos(x^4),  typeof(f(x)): Sym
  g:    sin(x^4) - cos(x^4),  typeof(g):    Sym

複数の数式は plot(), plot!() で図示する（一回の plot() 呼び出しだけで図示する方法はできなくなった）。

  f(x) = sqrt(x)
  fp = diff(f(x)) # fp は関数ではなく数式
  c = 2
  m = fp(x => c)  # at c=2
  println("c = $c,  m = $m")

  c = 2,  m = sqrt(2)/4

  plot(f(x), 1, 3)
  plot!(f(c) + m * (x - c))

5.1.　課題

5.1.1.　問題 1.

x^4 - 3x^3 + 3x^2 - 3x +2 を [-1, 3] の範囲で描画せよ。実数解はいくつあるか。

  g = x^4 - 3x^3 + 3x^2 - 3x + 2
  plot(g, -1, 3)
  savefig("fig4.png")

  solve(g) |> println
  real_roots(g) |> println
  real_roots(g) |> length |> println

  Sym[1, 2, -I, I]
  Sym[1, 2]
  2

6.　極限

limit() 関数は，数式の極限を求める。

limit() 関数の第 1 引数に数式を指定する。第 2 引数にどの変数について計算するか，第 3 引数はその変数をどこへ近づけていくかを指定する。

  limit(sin(x)/x, x, 0) |> println

  1

  limit((1-cos(x))/x^2, x, 0) |> println

  1/2

無限大は小文字の o を ２ 個続けて表す。

  limit(sin(x)/x, x, oo) |> println

  0

関数の導関数は limit() 関数を使って数式計算できる（数値計算もできるが）。

  @syms h
  f(x) = x^10
  limit((f(x + h) - f(x))/h, h, 0) |> println

  10*x^9

diff() での結果とあっている。

  diff(f(x)) |> println

  10*x^9

同じ結果を与えるがもっと複雑な計算式を使う方法もある。

  central_difference(f, x, h) = ( f(x + h) - f(x - h) ) / (2h)
  a = limit(central_difference(f, x, h), h, 0) |> println

  10*x^9

更にもっと複雑な計算式もある。

  central_4th_difference(f, x, h) = (-f(x + 2h) + 8f(x + h) - 8f(x - h) + f(x - 2h))/(12h)
  limit(central_4th_difference(f, x, h), h, 0) |> println

  10*x^9

limit() はチェイン規則にしたがう。

  func(x) = sin(x)
  limit((f(func(x + h)) - f(func(x))) / h, h, 0) |> println # 10*sin(1)^9*cos(1)

  10*sin(x)^9*cos(x)

6.1.　課題

6.1.1.　問題 1.

(3 - √(x + 5)) / (x - 4) の x を 4 に近づけていくときの極限を求めよ。

  f(x) = (3 - sqrt(x + 5)) / (x - 4)
  a = limit(f(x), x, 4)
  println("a = $a, float(a) = $(float(a))")

  a = -1/6, float(a) = -0.16666666666666666

6.1.2.　問題 2.

(x^(1/4) - 1) / (x^(1/3) - 1) の x を 1 に近づけていくときの極限を求めよ。

  f(x) = (x^(1/4) - 1) / (x^(1/3) - 1)
  limit(f(x), x, 1) |> float |> println

  0.75

6.1.3.　問題 3.

記号変数 h を定義し，f(x) = sin(x) とする。

  @syms h
  f(x) = sin(x);

(f(x+h) - 2f(x) + f(x-h)) / h^2 の h を 0 に近づけていくときの極限は，何になるか。

  gg = (f(x+h) - 2f(x) + f(x-h)) / h^2
  limit(gg, h, 0) |> println
  (diff(f(x), x), diff(f(x), x, 2), diff(f(x), x, 3)) |> println

  -sin(x)
  (cos(x), -sin(x), -cos(x))

この式は limit() 関数で 2 次の導関数を求めるための式である。

  f(x) = 3x^3 + 2x^2 +3x + 4
  limit(gg, h, 0) |> println
  expand(diff(f(x), x, 2)) |> println

  -sin(x)
  18*x + 4

6.1.4.　問題 4.

記号変数 a を定義し，

(cos(a * x) - 1) / (cos(x) - 1) において ｘ を 0 に近づけていくときの極限は，何になるか。

  @syms a
  limit((cos(a * x) - 1) / (cos(x) - 1), x, 0) |> println

  a^2

6.1.5.　問題 5.

式 f(x) = 1/(x^(log(log(log(1/x)))-1)) の x を 0 に近づけていくときの極限を数値的に求めることができる。

  f(x) = 1/(x^(log(log(log(1/x)))-1))
  xs = 0.1 .^ (2:4:40)
  fxs = [f(x) for x in xs]
  [xs fxs]

  10×2 Matrix{Float64}:
   0.01        0.0702822
   1.0e-6      0.619862
   1.0e-10    27.0054
   1.0e-14  2695.41
   1.0e-18     4.65934e5
   1.0e-22     1.20915e8
   1.0e-26     4.32215e10
   1.0e-30     2.00968e13
   1.0e-34     1.16712e16
   1.0e-38     8.21329e18

数式的に極限を求めると，結果は ∞ になり，数値的な場合に予測されたものに一致する（途中に float() を挟まない場合は，Inf を表す "oo" が表示される）。

  limit(f(x), x, 0) |> float |> println

  Inf

6.1.6.　問題 6.

式 f(x) = log(log(xexp(xexp(x)) + 1)) - exp(exp(log(log(x)) + 1/x)) で，x を ∞ に近づけていくときの極限を求めよ。

極限は 0 であるが，数値的に確認することはできない。

例えば f(10) をもとめてみよ。

  f(x) = log(log(x*exp(x*exp(x)) + 1)) - exp(exp(log(log(x)) + 1/x))
  f(x) |> println

  -exp(exp(1/x)*log(x)) + log(log(x*exp(x*exp(x)) + 1))

  f(10) |> println

  Inf

確かに x = 10 の時点で Inf になってしまっている。

しかし，x として長精度実数を使うと実は結果は Inf ではないことがわかる。

  f(big(10))

  -0.4374481647059431342589507553638095674342208743587679099109922069614458809810611

喜びもつかの間，x が大きくなるとそのうち（いくら使用する精度を高くしても） Inf になってしまう（要するに想像を絶する大きな数値になるのだ）。

  f(big(39))

  Inf

limit() 関数を使えば，正しい答え 0 が導かれる。

  limit(f(x), x, oo) |> println

  0

7.　導関数

前述したように，limit() を使って導関数を求めることもできるが，SymPy には導関数を求めるための diff() 関数がある。

以下で，f(x) は 関数ではなく，数式 Sym である。

  f(x) = exp(exp(x))

  f (generic function with 1 method)

数式中に変数が 2 個以上ある場合は，どの変数について導関数を求めるかを指定する（デフォルトは x）。

  diff(f(x)) |> println
  diff(f(x), x)  |> println

  exp(x)*exp(exp(x))
  exp(x)*exp(exp(x))

第三引数で何階導関数かを指定する（デフォルトでは 1）

  diff(f(x), x, 2) |> println

  (exp(x) + 1)*exp(x)*exp(exp(x))

x^x と，その導関数のグラフを描いてみる。

  f(x) = x^x
  plot(f(x), 1/10, 2, label="f(x)")
  plot!(diff(f(x)), label="f'(x)")

f(x) と c = 1 での切線のグラフを描く。

  f(x) = x^x; c = 1
  m = diff(f(x))
  plot(f(x), .5, 1.5)
  plot!(f(c) + m * (x - c))

7.1.　番外

diff() と solve() を組み合わせて，閉区間内の極値 f'(x) = 0 となる x を求めることができる。

3 辺の長さの和 L が一定の長方形の面積 A を最大にする問題を考えよう。

L = 2x + y
A = x * y

すなわち，A(x) = x * (L - 2x)

L を記号変数として，A'(x) = 0 となる x を求め，A(x) を求める。

  @syms L, x, y
  A = x*y
  A = A(y => L - 2x)
  out = solve(diff(A, x), x)[1]
  out |> println # x

  L/4

  (L - 2x)(x => out) |> println # y

  L/2

7.2.　課題

7.2.1.　問題 1.

tan(x)^(-1) の導関数を求めよ。

落とし穴（?）。tan(x)^(-1) は 1/tan(x) ではない。tan(x) の逆関数なので，atan(x) だ。

  diff(atan(x)) |> println

  1/(x^2 + 1)

微分は積分の逆なのだから，結果として提示されている 3 つの選択肢を integrate() してみればよい。atan(x) が解になる。

  integrate(1/(x^2 + 1)) |> println

  atan(x)

7.2.2.　問題 2.

f(x) = exp(-cos(x)) のとき，f'(3) を求めよ。

  diff(exp(-cos(x)))(x => 3) |> println 
  diff(exp(-cos(x)))(x => 3) |> float   |> println

  exp(-cos(3))*sin(3)
  0.3797841807457766324319619290395706689950768584935908050136628093384225996266999

7.2.3.　問題 3.

f(x) = 1/x^2 + x^3 の極値を求めよ。nsolve(diff(f(x)), 1) も試してみよ。

注: Julia では，? nsolve としても No documentation found. となるばかりだ。
Python の sympy で help(nsolve) でやっと，
Solve a nonlinear equation system numerically.
nsolve(f, [args,] x0, modules=['mpmath'])
'f' is a vector function of symbolic expressions representing the system.
args are the variables. If there is only one variable, this argument can
be omitted. 'x0' is a starting vector close to a solution.

  @syms x::real
  f(x) = 1/x^2 + x^3
  solve(diff(f(x)))              |> println
  solve(diff(f(x)))     |> float |> println
  nsolve(diff(f(x)), 1) |> float |> println # 1 は初期値

  Sym[2^(1/5)*3^(4/5)/3]
  [0.9221079114817278]
  0.9221079114817277656701705309291930146705850859540336029625316749125970624956978

7.2.4.　問題 4.

f(x) = tan(x) として， ニュートン法により 切線の方程式が f(c) + f'(c)(x-c) = 0 になるときの x を求めよ。

  f(x) = tan(x)
  @syms c
  solve(f(c) + diff(f(c)) * (x - c), x) |> println

  Sym[c - sin(2*c)/2]

c = pi/4 とすると何が得られるか？

  solve(f(c) + diff(f(c)) * (x - c), x)[1](c => PI/4) |> float |> println

  0.2853981633974483096156608458198757210492923498437764552437361480769541015715538

8.　積分（不定積分）

SymPy でも数式積分を行うことができるが，Mathematica に比べると劣る。

数値積分を行うのは integrate() である。

  @syms x, a
  f(x) = cos(x) - sin(x)
  integrate(f(x), x) |> println

  sin(x) + cos(x)

定数を含むこともできる。

  integrate(1/cos(x + a), x) |> println

  -log(tan(a/2 + x/2) - 1) + log(tan(a/2 + x/2) + 1)

その他の例

  integrate(x^3 * (1-x)^4, x) |> println

  x^8/8 - 4*x^7/7 + x^6 - 4*x^5/5 + x^4/4

  integrate(x/(x^2+1), x) |> println

  log(x^2 + 1)/2

9.　定積分

積分範囲を指定することで，定積分を計算することができる。

  integrate(x^2, (x, 0, 1)) |> println

  1/3

定積分の結果から，（分布関数の）中央値を計算することができる。

  using SymPy
  @vars x b real=true positive=true # 「b は正の実数」という条件を付けないと虚数解を含めて 4 個の解が求まる
  f(x) = 4x^3
  eq = integrate(f(x), (x,  0, b)) - 1/2 * integrate(f(x), (x, 0, 1));
  println(eq)
  solve(eq, b) |> float |> println

  b^4 - 0.5
  [0.8408964152537145]

9.1.　課題　　　

9.1.1.　問題 1.

f(x) = 2 / (x * log(x)) の不定積分を求めよ。

  integrate(2/(x * log(x))) |> println

  2*log(log(x))

9.1.2.　問題 2.

f(x) = (2x + 5)(x^2 + 5x)^7 を [0, 1] で定積分せよ。

  integrate((2x + 5) * (x^2 + 5x)^7, (x, 0, 1)) |> float |> println

  209952.0

9.1.3.　問題 3.

f(x) = sqrt(4 - sqrt(x)) を [0,9] で定積分せよ。

  integrate(sqrt(4 - sqrt(x)), (x, 0, 9)).evalf() |> println

  12.5333333333333

integrate(sqrt(4 - sqrt(x)), (x, 0, 9)) |> float |> println ではエラーになる

SymPy/Julia(1)

SymPy/Julia

https://www.math.csi.cuny.edu/~verzani/tmp/julia/symbolic.html

1.　Julia で数式処理

コンピュータによる計算処理で，３ 通りくらいのアプローチがある。

• Mathenatica，Maple, Sage に代表される数式処理
• MATLAB, R による数値処理
• Python, Haskel による一般化された汎用処理

Julia は 数値計算処理のために，MATLAB や R の代替として存在する。Julia は PyCall を使うことによって，Python およびその SymPｙ パッケージへの連携を提供した。

Julia で SymPy パッケージを使えるようにする方法については，ここでは述べない。単に，Python を動くようにすることと，SymPy パッケージを使えるようにした後，Pkg.add("SymPy") するだけである。

2.　概略

数式処理をするために必用なのは，SymPy を使うことを宣言する（using SymPy）と，そこで使う変数を宣言すること（@syms x など）である。

  using SymPy
  @syms x

  (x,)

これによって宣言（定義）される ｘ は，Julia における変数名ではなく，以下に引き続く処理の記述で使用される記号変数（Sym オブジェクト）である。

x の型(typeof) は Sym である。

  x |> typeof |> println

  Sym

Sym オブジェクト（以下では単に変数と呼ぶ）は，何通りかの方法で宣言（定義）することができるが，@syms マクロ（または @vars マクロ）を使うのがもっと便利であろう。h, y の 2 変数を定義するには以下のようにする。

変数名は空白あるいはカンマで区切って複数記述できる。

  @syms h y

  (h, y)

なお，実際の使用例を後で示すが，特定の条件をもつ変数として宣言することもできる。

  @syms a::positive b::real c::integer

  (a, b, c)

複数の条件を持つ場合は @vars を使い，条件をタプルで指定する。

  @vars x y z real=true w nonnegative=true v (positive=true, real=true)

  (x, y, z, w, v)

このようにして定義した変数を使用して，任意の数式を記述することができる。

たとえば x^2 -2x +2 のような数式はそのまま指定できる。

この記事を通じて，数式のあとに "|> typeof" あるいは "|> println" としているのは，左側の数式に対して typeof(数式)，さらに printiln(処理された数式） を実行し，最終的な結果を表示する（REPL の場合）ことを示しているだけである。

ユーザが実際に SymPy を使用する環境によってはそれらは不要，あるいはないほうがよい，更にはあっては困るということもあるであろうから，適宜省略・削除するようお願いしたい。（最後のいくつかの |> を省略したらどうなるか，事前に把握しておいて欲しい）。

以下では，定義した数式の型が何であるかを表示している。なれてくれば，このような確認は不用であろう。

  x^2 - 2x + 2 |> typeof |> println

  Sym

以下では，関数名 f で普通の Julia 関数を定義している。

  f(x) = exp(-x^2/2)

  f (generic function with 1 method)

普通通り，x = 3 のときの関数値は f(3) で得られる。

  f(3) |> println

  0.011108996538242306

しかし，f(x) は Sym 型を持つ数式である（typeof(f(x)) = Sym）。

  f(x) |> println

  exp(-x^2/2)

  typeof(f(x)) |> println

  Sym

f(x) が Sym オブジェクトであるということは，そのオブジェクト中の変数を特定の値として評価できるということである。

以下は，f(x) 中の x を 1 にするとどのような数式になるかを示すものである。

  subs(f(x), x, 1) |> println

  exp(-1/2)

  f(x)(x => 1) |> println

  exp(-1/2)

f(x) が表す式にある x を 1 で置き換えると exp(-1/2) になるということである。

subs() 関数によって式の一部が置き換えられた（代入された）ものも，前と同じ Sym オブジェクトなので，その式に別のものを代入することもできる。

以下の示すのは，f(x) において， x = 1 を代入 (subs) した結果の数値 (float) を表示する (println) という一連の処理を示している。

  f(x) |> subs(x => 1) |> float |> println

  0.6065306597126334

同じ処理をするのには何通りかの方法があることがあるが，今の場合は以下のほうがより簡潔かもしれなお。

  f(x)(x => 1) |> float |> println

  0.6065306597126334

多くの場合，最終的な（人間が見てわかる）表示形式にするのは一番最後になるであろう。

3.　方程式

SymPy は代数演算として，簡約化，因数分解，方程式の解など多くのタスクをこなすことができる。

しかし，数式処理をするために必用とされるテクニックはそんなに多くはない。

例えば，以下のように記述すれば，数式は自動的に簡約化される（最後の 「|> println」 は，処理された結果を画面表示するとそのようになるということであるが，環境によってはそれがないほうがわかりやすい数式表現がされることもあろう）。

  (x+1) + (x+2) + (x+3) |> println

  3*x + 6

書かれた数式をもう少しわかりやすい（人と場合によるのだが）形式で表現してほしいこともある。

例えば，以下の式は人によっては別の形式で表現してほしいかもしれない。現役数学者（中高生？）なら，我慢できないかもしれない。

  sin(x)^2 + cos(x)^2 |> println

  sin(x)^2 + cos(x)^2

このよう場合，simplyfi() 関数により有名すぎる公式「x が何であろうと（有効範囲はあるが），sin(x) と cos(x) の二乗の和和は 1」簡約化された解を得られる。

  simplify(sin(x)^2 + cos(x)^2) |> println

  1

  sin(x)^2 + cos(x)^2 |> simplify |> println

  1

factor() 関数により，数式の因数分解もできる。

  factor(x^2 + 2x + 1) |> println

  (x + 1)^2

因数分解された形式で表示された数式は，expand() 関数によって展開された数式にすることができる。

  expand( (x+1)*(x+2)*(x+3) ) |> println

  x^3 + 6*x^2 + 11*x + 6

4.　方程式を解く

多項式の実数解を数式解として求めるには real_roots() 関数を使う。

  real_roots(x^2 + x - 1) |> println

  Sym[-sqrt(5)/2 - 1/2, -1/2 + sqrt(5)/2]

もう少し一般的にいえば，多項式に限らず方程式の数式解を得るには solve() 関数を使う。

  solve(x^2 + x -1, x) |> println

  Sym[-1/2 + sqrt(5)/2]

得られるのは数式解なので，数値解を得るには float() 関数を適用すればよい。

  solve(x^2 + x -1, x) |> float |> println

  [0.6180339887498949]

三角関数方程式は無数の解を持つが，SymPy の sin では -90度 〜 90度の範囲の解を求める。

  solve(sin(x)^2 - (1/2)^2) * (180 /pi) |> println

  Sym[30.0000000000000, 150.000000000000, 210.000000000000]

式中に変数が 2 つ以上ある場合に，どの変数についての解を求めるかは第 2 引数で指定する。

  fx = x^2 + y^2 - 1
  out = solve(fx, y)
  println(out)

  Sym[-sqrt(1 - x^2), sqrt(1 - x^2)]

方程式は二次式なので，2 つの数式解を含むベクトルが返される。x = 1/2 のときの値は，各要素について x = 1/2 を代入して，float にする。

  out[1](x => 1/2) |> float |> println
  out[2](x => 1/2) |> float |> println

  -0.8660254037844386
  0.8660254037844386

以下のようにすればベクトルとしてまとめて表示できる。

  println([o(x => 1/2) |> float for o in out])

  [-0.8660254037844386, 0.8660254037844386]

  println(map(arg -> arg(x => 1/2) |> float, out))

  [-0.8660254037844386, 0.8660254037844386]

連立方程式も solve で解くことができる。 x + y - 2a = 0 x - y + b = 0 を x, y について解く。 連立方程式と，どの変数について解くかの指定をベクトルで指定する。

  solve([x + y - 2a, x - y + b], [x, y]) |> println

  Dict{Any, Any}(y => a + b/2, x => a - b/2)

別法として，Julia で数値解を得るときと同じように記述することもできる。

x + y = 2a x - y = -b

係数行列 [1 1; 1 -1] と定数ベクトル [2a, -b] と \ 演算子を使う。

  Sym[1 1; 1 -1] \ Sym[2a, -b] |> println

  Sym[a - b/2, a + b/2]

4.1.　課題

4.1.1.　問題 1.

(x-1)(x-2)(x-3)(x-4)(x-5) において，x^3 の係数はいくつか。

  eq = expand((x-1)*(x-2)*(x-3)*(x-4)*(x-5))
  eq |> println

  x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120

目で見るだけで解は分かるが，係数を取り出す object.coeff(...) メソッドを使えばよい。

  eq.coeff(x^3) |> println

  85

4.1.2.　問題 2. Horner 法

((1x + 2)x + 3)x + 4 を展開せよ。

  ((2x + 3)x + 4)x + 5 |> expand |> println

  2*x^3 + 3*x^2 + 4*x + 5

Horner 法は，多項式 2 x^3 + 3 x^2 + 4 x + 5 を記述する方法の一つで，((2 x + 3)x + 4)x + 5 のように ( ) を使って記述する。

((2 x + 3) x + 4) x + 5 では，加算が 3 回，乗算が 3 回であるのに対し
2 x^3 + 3 x^2 + 4 x + 5 では，加算は同じく 3 回であるが，乗算は（^ も加え） 6 回（3, 2, 1 回）である。 もっとも，最近のコンピュータの性能から見ればたいした違いではないだろうし，Julia のようなコンパイラ言語では数式の解釈時に Hoener 法を使って最適化するようなこともあるのかもしれない。

4.1.3.　問題 3.

x^3 - 6x^2 + 11x - 6 = 0 を解き，最も大きい解を求めよ。

以下で，N() 関数は，適切な数値表現への変換関数である。float() 関数は引数が整数でも実数化してしまう（この場合には，実害はないのだが）。

  solve(x^3 - 6x^2 + 11x - 6) |> N |> println

  [1, 2, 3]

目で見ても分かるが，最大値を求める関数 maximum() を使えばよい。

  solve(x^3 - 6x^2 + 11x - 6) |> maximum |> println

  3

4.1.4.　問題 4.

直角三角形において，直角を挟む二辺が 3cm と 4cm の場合に，斜辺の長さ x を求めよ。

この場合，x は非負の実数なので，@syms で変数を定義するときに「変数::変数の特性」の書き方を使い，nonnnegative であるという条件を付けておくと不要な解が含まれなくなる。

  @syms x::nonnegative
  solve(x^2 - (3^2 + 4^2)) |> N |> println

  [5]

なお，Julia ではベクトルは要素が １ 個の場合でもベクトルなので，スカラーとして取り出すときは 「ベクトル[1]」 のように添え字を指定しないと「[要素]」となってしまうので注意が必要である。

要素数が 1 個の場合と，2 個以上の場合の両方に対応するには以下のようにする。

  a = solve(x^2 - (3^2 + 4^2))
  length(a) == 1 ? println(N(a[1])) : println(a)

  5

  @syms y
  b = solve((y + 3) * (y - 5))
  length(b) == 1 ? println(N(b[1])) : println(N(b))

  [-3, 5]

4.1.5.　問題 5.

x^4 - 10x^3 + 32x^2 - 38x + 15 = 0 を解き，最も大きい解を求めよ。

問題 5 と同じである。

  @syms x
  solve(x^4 - 10x^3 + 32x^2 - 38x + 15) |> N |> println

  [1, 3, 5]

  solve(x^4 - 10x^3 + 32x^2 - 38x + 15) |> maximum |> println

  5

いまの場合，4 次式なので解は 4 個あるはずだが，solve() では 3 個しか表示されなかった。

solve では重解は複数として数えないのである。factor を使えば詳細が分かる。

  factor(x^4 - 10x^3 + 32x^2 - 38x + 15) |> println

  (x - 5)*(x - 3)*(x - 1)^2

この場合は，x = 1 が二重解であることが分かる。

このような場合には，roots() 関数を使えば重解も複数として数える。解は辞書形式で返される。

  roots(x^4 - 10x^3 + 32x^2 - 38x + 15) |> println

  Dict{Any, Any}(3 => 1, 1 => 2, 5 => 1)

この場合は，x=3, x=5 は 1つずつ，x=1 は 2 個であることが分かる。

実数解だけを求める必要がある場合は，real_roots() 関数を使うことができる。

以下の p は 8 次式なので，実数解，虚数解を合わせて 8 個ある。

  p = (x - 3)^2 * (x - 2) * (x-1) * x * (x + 1) * (x^2 + x + 1)
  expand(p) |> println

  x^8 - 7*x^7 + 13*x^6 + 2*x^5 - 11*x^4 - 13*x^3 - 3*x^2 + 18*x

real_roots() によれば，そのうち 6 個は実数解であることがわかる。

  real_roots(p) |> N |> println

  [-1, 0, 1, 2, 3, 3]

p のすべての解は roots() 関数でわかる。

虚数解は共役複素数 -1/2 ± sqrt(3)*I/2 （2 個）である。

  roots(p) |> println

  Dict{Any, Any}(-1 => 1, 3 => 2, 1 => 1, -1/2 - sqrt(3)*I/2 => 1, -1/2 + sqrt(3)*I/2 => 1, 0 => 1, 2 => 1)

4.1.6.　問題 6.

e^x = x^4 を解け。

Julia で，ℯ は REPL で \euler[tab] ，⩵ は \Equal[tab] として入力する。（ここからコピーペーストしても良い）

  solve(ℯ^x ⩵ x^4) |> println

  Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]

⩵ が入力しづらければ Eq() 関数を使って以下のように書くこともできる。

  solve(Eq(ℯ^x, x^4)) |> println

  Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]

LambertW() Lambert の W 関数というものだそうだ。float() を通せばその値がわかる。

  solve(Eq(ℯ^x, x^4)) |> float |> println

  [1.4296118247255556, -0.8155534188089607, 8.6131694564414]

N() では BigFloat として多くの桁数が表示される

  solve(Eq(ℯ^x, x^4)) |> N |> println

  BigFloat[1.42961182472555561227524441622361901326636222030404817451048179435868561038458, -0.8155534188089606577727273253085594805974099088406385398936250403272774562683745, 8.613169456441398596676396600371925502144825941277910783538046030885450338602536]

方程式の解として正しいか，evalf() を使って検算してみる。

3 番目の解では，左辺と右辺の数値は表示されている範囲では等しいが，全体としては誤差の範囲で一致しないということであろう。

  for res in solve(Eq(ℯ^x, x^4))
      left, right = ℯ^res.evalf(), (res^4).evalf()
      println("left = $left  right = $right  right == left = $(right == left)")
  end

  left = 4.17707743900016  right = 4.17707743900016  right == left = true
  left = 0.442394430203601  right = 0.442394430203601  right == left = true
  left = 5503.66468907673  right = 5503.66468907672  right == left = false

Lambert の W 関数は Julia の LambertW パッケージに用意されている。

  # import Pkg; Pkg.add("LambertW")
  using LambertW
  solve(Eq(ℯ^x, x^4)) |> println

  Sym[-4*LambertW(-1/4), -4*LambertW(1/4), -4*LambertW(-1/4, -1)]

  typeof(sympy.LambertW)

  PyObject

PyCall.PyObject の LambertW() 関数を LambertW パッケージの lambertw() 関数で置き換える。

  res = [-4*lambertw(-1/4), -4*lambertw(1/4), -4*lambertw(-1/4, -1)]
  res |> float |> println

  [1.4296118247255556, -0.8155534188089607, 8.6131694564414]

検算してみる。

  res .^ 4 |> float |> println

  [4.177077439000158, 0.44239443020360086, 5503.664689076726]

  ℯ .^ res |> float |> println

  [4.177077439000158, 0.44239443020360075, 5503.664689076728]