Edit on GitHub

第9回： ■ 配列要素の操作／▶常微分方程式の数値解法

■ ベクトルを引数とする関数

■ 総和関数 sum のように，ベクトルを引数とする関数がある．

■ 積

julia> v = [2, 3, 4];

julia> prod(v)
24

julia> r = 1;

julia> for i = 1:length(v)
          global r
          r *= v[i]
       end

julia> r
24

■ ノルム

「ノルム（norm）」は，ベクトル（や行列）の「大きさ」を一般化した関数である．

LinearAlgebra パッケージの中で，関数 norm() が定義されている．

ノルムにはいくつかの定義がある． 単なる norm(v) は，2-norm を意味し，各要素の２乗平均値の和の平方根である．

julia> v = [1, 2, 3, 4, 5, 6, 7];

julia> using LinearAlgebra

julia> norm(v)
11.832159566199232

julia> @show sqrt(sum(v .^ 2))
sqrt(sum(v .^ 2)) = 11.832159566199232
11.832159566199232

julia> r = 0;

julia> for i = 1:length(v)
          global r
          r += v[i]^2
       end

julia> sqrt(r)
11.832159566199232
Note

関数 abs.(v) は，ベクトルの各要素の絶対値からなるベクトルである．

julia> v = [1, 2, 3, 4, 5, 6, 7];

julia> abs.(v)
7-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6
 7

■ 平均値・標準偏差

ベクトルに格納されたデータの平均値や標準偏差を計算できる．

Statistics パッケージの関数 mean(v) は，ベクトル v の平均値を算出する．平均値は，各要素の総和 sum(v) を要素の数 $n$ で除したものである．

Statistics パッケージの関数 std(v) は，ベクトル v の標準偏差を算出する．

単なる std(v) は，$(n-1)$ で割った「偏りがない（unbiased）」標準偏差を算出する． 平均値を算出する． std(v, corrected=false) とすると，$n$ で割った「偏った（biased）」標準偏差を算出する．

julia> v = [1, 2, 3, 4, 5, 6, 7];

julia> # Statistics パッケージの読み込み
       using Statistics

julia> # 平均値
       mean(v)
4.0

julia> sum(v) / length(v)
4.0

julia> # 偏りがない標準分散，n-1 で割る
       std(v)
2.160246899469287

julia> sqrt(sum((v .- mean(v)) .^ 2) / (length(v) - 1))
2.160246899469287

julia> # 偏った標準分散，n で割る
       std(v, corrected = false)
2.0

julia> sqrt(sum((v .- mean(v)) .^ 2) / (length(v)))
2.0
Note

標準分散の計算には，「偏りのない」定義を用いるのがよい．例えば，こちらを参照．→ 分散は n で割るか n − 1 で割るか

■ 複数の数を引数とする関数

julia> min(5, 1, 4, 2, 3)
1

julia> max(5, 1, 4, 2, 3)
5

■ splatting 演算子

...演算子は，関数呼び出しにおいて，ベクトルを，複数の引数に分けてから呼び出す．

julia> min([5, 1, 4, 2, 3]) # => exception
ERROR: MethodError: no method matching min(::Array{Int64,1})
Closest candidates are:
  min(::Any, !Matched::Missing) at missing.jl:126
  min(::Any, !Matched::Any) at operators.jl:431
  min(::Any, !Matched::Any, !Matched::Any, !Matched::Any...) at operators.jl:538
  ...

julia> min([5, 1, 4, 2, 3]...) # min(5,1,4,2,3) と同じ
1

■ ベクトル要素への代入

julia> v = collect(1:10)
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> # インデックス：整数
       v[4] = 0
0

julia> v
10-element Array{Int64,1}:
  1
  2
  3
  0
  5
  6
  7
  8
  9
 10

演算子 .= は，ベクトルの各要素に対する代入である． ベクトルの要素を，整数の等差数列で指定して，一度に更新できる．

julia> # インデックス：範囲
       v[3:2:10] .= 0
4-element view(::Array{Int64,1}, 3:2:9) with eltype Int64:
 0
 0
 0
 0

julia> v
10-element Array{Int64,1}:
  1
  2
  0
  0
  0
  6
  0
  8
  0
 10

julia> # `=` では例外を発生する
       v[3:2:10] = 1 # => Exception
ERROR: ArgumentError: indexed assignment with a single value to many locations is not supported; perhaps use broadcasting `.=` instead?

■ 素数の生成：エラトステネスの篩

エラトステネスの篩（ふるい）は，素数を算出する方法の一つである． 以下の手順による．

  • 数 $2$ から $n$ までの整数を並べる
  • 生き残っている中で最も小さい数 $p$ を素数として残す．
  • 素数 $p$ 自身を除く $p$ の倍数を全て消す
  • 以上の手順を，$n$ まで調べたら終わり．

以下のプログラムでは，配列 sieve を篩とする． 篩の初期値を 1:n とすると，数字 i の篩は sieve[i] である． 篩で消された数 $i$ には sieve[i] に 0 を格納することにする．

nmax = 100
sieve = collect(1:nmax);
sieve[1] = 0;
for i = 2:nmax
   if sieve[i] > 0
      println(i)
      for j = i*2:i:nmax
         sieve[j] = 0
      end
   end
end
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97

上のプログラムで，変数 j に関する繰り返しは，1行で書ける．

nmax = 100
sieve = collect(1:nmax);
sieve[1] = 0;
for i = 2:nmax
   if sieve[i] > 0
      # println(i)
      sieve[i*2:i:nmax] .= 0
   end
end

for i = 1:nmax
   if sieve[i] > 0
      println(i)
   end
end
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97

ここで，sieve[i*2:i:nmax].=0 の文は， 等差数列 i*2:i:nmax で表される添字で示される配列 sieve の全てに 0 を代入することを意味する． 等差数列 i*2:i:nmax は，i*2から始まり，i 飛びに nmax まで増える等差数列である．

Note

Julia には，素数を高速に計算する関数を含むパッケージが用意されている．

primes(n) は，数 $n$ までの素数を計算する．

isprime(x)は，数 $x$ が素数であるかどうかを判定する．

julia> # import Pkg; Pkg.add("Primes") # コメントを外してパッケージを設置する．一度だけ行えばよい
       using Primes

julia> isprime(2)
true

julia> isprime(3)
true

julia> isprime(4)
false

julia> isprime.([2, 3, 4])
3-element BitArray{1}:
 1
 1
 0

julia> primes(100)
25-element Array{Int64,1}:
  2
  3
  5
  7
 11
 13
 17
 19
 23
 29
  ⋮
 59
 61
 67
 71
 73
 79
 83
 89
 97

▶ 常微分方程式の初期値問題

▶ 常微分方程式の初期値問題:Euler法

微分方程式

\[\dfrac{dx}{dt} =f(x,t),\]

の解 $x(t)$ （の近似値）を求めたい．

Euler 法による数値解法は，以下のような手順である．

時刻 $t_1, t_2, \ldots$ を一定間隔 $h$ とする． 上の式を，以下のように離散化する．

\[\begin{aligned} \dfrac{x_{n+1}-x_{n}}{h} & = f(x_n,t_n) \\ x_{n+1} & = x_{n} + h f(x_n,t_n) \end{aligned}\]

例題：Euler法による解法

以下の微分方程式を解いてみる．

\[\begin{aligned} \dfrac{dx}{dt} & = 1-x^2, \\ x(0) & = 0, \\ 0 & \leq t \leq 1.6 \end{aligned}\]

刻み $h = 0.4$ とする．

f(x, t) = 1 - x^2
#
tmin = 0
tmax = 1.6
h = 0.4
ts = tmin:h:tmax
n = length(ts)
#
x_now = 0 # initial condition
for i = 1:n
   global x_now
   t = ts[i]
   x_next = x_now + h * f(x_now, t)
   @show t, x_next
   x_now = x_next
end
(t, x_next) = (0.0, 0.4)
(t, x_next) = (0.4, 0.736)
(t, x_next) = (0.8, 0.9193216)
(t, x_next) = (1.2, 0.981260718309376)
(t, x_next) = (1.6, 0.996111679390563)

解析解は，$x = \tanh{t}$ である．

using PyPlot

x_now = 0 # initial condition
for i = 1:n
   global x_now
   t = ts[i]
   plt.plot(t, x_now, "b.")
   x_next = x_now + h * f(x_now, t)
   @show t, x_next
   x_now = x_next
end
plt.plot(ts, tanh.(ts), "r")
plt.xlabel("t")
plt.ylabel("x")
(t, x_next) = (0.0, 0.4)
(t, x_next) = (0.4, 0.736)
(t, x_next) = (0.8, 0.9193216)
(t, x_next) = (1.2, 0.981260718309376)
(t, x_next) = (1.6, 0.996111679390563)

配列に計算結果を入れて，一気に描画する．

using PyPlot

tmin = 0
tmax = 1.6
h = 0.4
ts = tmin:h:tmax
n = length(ts)
xs = zeros(n)
xs[1] = 0 # initial condition

for i = 1:n-1
   local x_now = xs[i]
   t = ts[i]
   x_next = x_now + h * f(x_now, t)
   xs[i+1] = x_next
end
plt.plot(ts, xs, ".")
plt.plot(ts, tanh.(ts), "r")
plt.xlabel("t")
plt.ylabel("x")

刻みを狭くする

刻み $h$ を $0.4, 0.2, 0.1, 0.05$ と小さくしてみる． 刻みを小さくすると，近似解が厳密解に近づいていくことが観察できる．

using PyPlot

tmin = 0
tmax = 1.6
h = 0.4
for k = 1:4
   global h
   local ts = tmin:h:tmax
   local n = length(ts)
   local xs = zeros(n)
   xs[1] = 0 #  initial condition

   for i = 1:n-1
      t = ts[i]
      local x_now = xs[i]
      x_next = x_now + h * f(x_now, t)
      xs[i+1] = x_next
   end
   plt.plot(ts, xs, ".", label = "h=" * string(h))

   h /= 2
end
plt.xlabel("t")
plt.ylabel("x")
plt.plot(ts, tanh.(ts), "b", label = "tanh(t)", lw = 0.5)
plt.legend()

正確な解との誤差評価

using LinearAlgebra
using PyPlot

tmin = 0
tmax = 1.6
h = 0.4
for k = 1:5
   global h
   local ts = tmin:h:tmax
   local n = length(ts)
   local xs = zeros(n)
   xs[1] = 0 #  initial condition

   for i = 1:n-1
      t = ts[i]
      local x_now = xs[i]
      x_next = x_now + h * f(x_now, t)
      xs[i+1] = x_next
   end
   xtrue = tanh.(ts)
   e = norm(xs .- xtrue) / n
   @show h, e
   plt.plot(h, e, ".")
   h /= 2
end
plt.xlabel("h")
plt.xscale("log")
plt.yscale("log")
plt.xlim(1e-2, 1)
plt.ylim(1e-4, 1e-1)
(h, e) = (0.4, 0.025667730896465946)
(h, e) = (0.2, 0.00931239766406314)
(h, e) = (0.1, 0.0033516722128539393)
(h, e) = (0.05, 0.0011971170296258557)
(h, e) = (0.025, 0.00042554173573778107)

▶ 常微分方程式の初期値問題:修正Euler法

修正Euler 法では，微分方程式

\[\dfrac{dx}{dt} =f(x,t)\]

を，次のように離散化する．

\[\begin{aligned} f_{n} & = f(x_{n}, t_{n}), \\ \overline{x}_{n+1} & = x_{n} + h f(x_n,t), \\ \overline{f}_{n+1} & = f(\overline{x}_{n+1}, t_{n+1}) \\ x_{n+1} & = x_{n} + \dfrac{h}{2} \left(f_{n} + \overline{f}_{n+1}\right) \end{aligned}\]

例題：修正Euler法による解法

（再掲）Euler法と同じ微分方程式を解いてみる．

\[\begin{aligned} \dfrac{dx}{dt} & = 1-x^2, \\ x(0) & = 0, \\ 0 & \leq t \leq 1.6 \end{aligned}\]

刻み $h = 0.4$ とする．

#
tmin = 0
tmax = 1.6
h = 0.4
ts = tmin:h:tmax
n = length(ts)
x_now = 0 # initial condition

for i = 1:n-1
   global x_now
   t = ts[i]
   t_next = ts[i+1]
   f_now = f(x_now, t)
   x_mid = x_now + h * f_now
   f_mid = f(x_mid, t_next)
   x_next = x_now + (f_now + f_mid) * h / 2
   @show t, x_next
   x_now = x_next
end
(t, x_next) = (0.0, 0.368)
(t, x_next) = (0.4, 0.6390044320071679)
(t, x_next) = (0.8, 0.8039781901649501)
(t, x_next) = (1.2, 0.8959360086216626)

配列に計算結果を入れて，一気に描画する．

using PyPlot

n = length(ts)
xs = zeros(n)

xs[1] = 0 # initial condition
for i = 1:n-1
   global xs
   local x_now = xs[i]
   t_next = ts[i+1]
   t = ts[i]
   f_now = f(x_now, t)
   x_mid = x_now + h * f_now
   f_mid = f(x_mid, t_next)
   xs[i+1] = x_now + (f_now + f_mid) * h / 2
end
plt.plot(ts, xs, ".")
plt.plot(ts, tanh.(ts))
plt.xlabel("t")
plt.ylabel("x")

刻みを狭くする

刻み $h$ を $0.4, 0.2, 0.1, 0.05$ と小さくしてみる． 刻みを小さくすると，近似解が厳密解に近づいていくことが観察できる．

using LinearAlgebra
using PyPlot

h = 0.4
for k = 1:4
   global h
   local ts = tmin:h:tmax
   local n = length(ts)
   local xs = zeros(n)
   xs[1] = 0 # initial condition
   for i = 1:n-1
      t = ts[i]
      local x_now = xs[i]
      t_next = ts[i+1]
      f_now = f(x_now, t)
      x_mid = x_now + h * f_now
      f_mid = f(x_mid, t_next)
      xs[i+1] = x_now + (f_now + f_mid) * h / 2
   end
   xtrue = tanh.(ts)
   e = norm(xs .- xtrue)
   @show h, e
   plt.plot(ts, xs, ".", label = "h=" * string(h))
   h /= 2
end
plt.xlabel("t")
plt.ylabel("x")
plt.plot(ts, tanh.(ts), "b", label = "tanh(t)", lw = 0.5)
plt.legend()
(h, e) = (0.4, 0.048085853296269084)
(h, e) = (0.2, 0.013351045559265527)
(h, e) = (0.1, 0.0042050468178388605)
(h, e) = (0.05, 0.001404770260316803)

正確な解との誤差評価

using LinearAlgebra
using PyPlot

h = 0.4
for k = 1:4
   global h
   local ts = tmin:h:tmax
   local n = length(ts)
   local xs = zeros(n)
   xs[1] = 0 # initial condition
   for i = 1:n-1
      t = ts[i]
      local x_now = xs[i]
      t_next = ts[i+1]
      f_now = f(x_now, t)
      x_mid = x_now + h * f_now
      f_mid = f(x_mid, t_next)
      xs[i+1] = x_now + (f_now + f_mid) * h / 2
   end
   xtrue = tanh.(ts)
   e = norm(xs .- xtrue) / n
   @show h, e
   plt.plot(h, e, ".")
   h /= 2
end
plt.xlabel("h")
plt.xscale("log")
plt.yscale("log")
plt.xlim(1e-2, 1)
plt.ylim(1e-5, 1e-1)
(h, e) = (0.4, 0.009617170659253816)
(h, e) = (0.2, 0.0014834495065850586)
(h, e) = (0.1, 0.0002473556951669918)
(h, e) = (0.05, 4.2568795767175844e-5)

◀● 練習：常微分方程式の数値解の誤差

上の常微分方程式の数値解法の例について， Euler法による絶対誤差と，修正Euler法による絶対誤差を， 刻み幅 $h$ に対する関数として，一つのグラフの上に表せ．

結果は，例えば，以下のようになろう．

◀● 練習： 条件が成り立つまで繰り返す：微分方程式の初期値問題

(少し難しいので，後回しにしてもよい)

Euler法ないし修正Euler法による微分方程式の数値解法を， 刻み幅 $h$ を半分にしながら 20 回繰り返せ． ただし，絶対誤差が $10^{-4}$ 以下になったら，そこで終了せよ．

◀● 練習：常微分方程式・素性の悪い問題

以下の微分方程式を解いてみよ．

\[\begin{aligned} \dfrac{dx}{dt} & = x^2, \\ x(0) & = \dfrac{1}{2}, \\ 0 & \le t < 2 \end{aligned}\]

解析解は，

\[x = \dfrac{1}{2-t}\]

となり，$t \longrightarrow 0$ で無限大に発散する「素性の悪い」方程式である．

★ 今回のまとめ

  • ベクトルを引数とする関数
  • 複数の数を引数とする関数
  • splatting 演算子
  • ベクトル要素への代入
  • エラトステネスの篩：素数を算出する
  • 微分方程式の初期値問題，Euler法，修正Euler法