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

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

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

■ 積

• Base.prod — Function

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> r24

■ ノルム

• LinearAlgebra.norm — Function

「ノルム（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

julia> v = [1, 2, 3, 4, 5, 6, 7];
julia> abs.(v)7-element Vector{Int64}:
 1
 2
 3
 4
 5
 6
 7

■ 平均値・標準偏差

• Statistics.mean — Function
• Statistics.stdm — Function

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

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

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

• Base.min — Function
• Base.max — Function

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

■ splatting 演算子

• ... splits one argument into many different arguments in function calls

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

julia> min([5, 1, 4, 2, 3]) # => exceptionERROR: MethodError: no method matching min(::Vector{Int64})
Closest candidates are:
  min(::Any, !Matched::Missing) at missing.jl:134
  min(::Any, !Matched::Any) at operators.jl:494
  min(::Any, !Matched::Any, !Matched::Any, !Matched::Any...) at operators.jl:591
  ...
julia> min([5, 1, 4, 2, 3]...) # min(5,1,4,2,3) と同じ1

■ ベクトル要素への代入

julia> v = collect(1:10)10-element Vector{Int64}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
julia> # インデックス：整数
       v[4] = 00
julia> v10-element Vector{Int64}:
  1
  2
  3
  0
  5
  6
  7
  8
  9
 10

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

julia> # インデックス：範囲
       v[3:2:10] .= 04-element view(::Vector{Int64}, 3:2:9) with eltype Int64:
 0
 0
 0
 0
julia> v10-element Vector{Int64}:
  1
  2
  0
  0
  0
  6
  0
  8
  0
 10
julia> # `=` では例外を発生する
       v[3:2:10] = 1 # => ExceptionERROR: ArgumentError: indexed assignment with a single value to possibly 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 まで増える等差数列である．

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 BitVector:
 1
 1
 0
julia> primes(100)25-element Vector{Int64}:
  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.025667730896465974)
(h, e) = (0.2, 0.00931239766406315)
(h, e) = (0.1, 0.0033516722128539445)
(h, e) = (0.05, 0.001197117029625857)
(h, e) = (0.025, 0.000425541735737782)

▶ 常微分方程式の初期値問題:修正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.048085853296268966)
(h, e) = (0.2, 0.013351045559265437)
(h, e) = (0.1, 0.004205046817838765)
(h, e) = (0.05, 0.0014047702603167624)

正確な解との誤差評価

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.009617170659253793)
(h, e) = (0.2, 0.0014834495065850486)
(h, e) = (0.1, 0.0002473556951669862)
(h, e) = (0.05, 4.256879576717462e-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法