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 ではエラーになる