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

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

■ 総和関数 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

Note

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

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(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];
       # Statistics パッケージの読み込み
julia> using Statistics
       # 平均値
julia> mean(v)4.0
julia> sum(v) / length(v)
       # 偏りがない標準分散，n-1 で割る4.0
julia> std(v)2.160246899469287
julia> sqrt(sum((v .- mean(v)) .^ 2) / (length(v) - 1))
       # 偏った標準分散，n で割る2.160246899469287
julia> 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` 演算子

... 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, ::Missing)
   @ Base missing.jl:134
  min(::Any, ::Any)
   @ Base operators.jl:490
  min(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:587
  ...
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> v
       # `=` では例外を発生する10-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?

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

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

数 $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

上のプログラムで，変数 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

ここで，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 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

プロットデフォルトの設定10

using LinearAlgebra

aux = (label="", lc=:black, lw=0.5)
using Plots
using Plots.PlotMeasures

Plots.reset_defaults()
default(
    size=(720,480),
    framestyle=:box,thickness_scaling=2,
    dpi=150,
    guidefontsize=10,
    tickfontsize=8,
    legendfontsize=8,
    aspect_ratio=:auto,
    label="",
    legend_foreground_color=nothing,
    legend_background_color=nothing,
    grid=false,xlabel="", ylabel="",
    palette=:Set1_9,
    linewidth=2,
    markershape=:none,
    markersize=2,
    markerstrokewidth=0,
    warn_on_unsupported=true)

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

▶ 常微分方程式の初期値問題: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 Plots
plot(xlabel="t", ylabel="x")

x_now = 0 # initial condition
for i = 1:n
   global x_now
   t = ts[i]
   scatter!([t], [x_now], mc=1, shape=:circle)
   x_next = x_now + h * f(x_now, t)
   @show t, x_next
   x_now = x_next
end
plot!(ts, tanh.(ts), lc=2)

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

plot(xlabel="t", ylabel="x")

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
scatter!(ts, xs, mc=1, shape=:circle)
plot!(ts, tanh.(ts), lc=2)

刻みを狭くする

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

plot(xlabel="t", ylabel="x")

tmin = 0
tmax = 1.6
h = 0.4
for k = 1:4
   global h
   ts = tmin:h:tmax
   n = length(ts)
   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
   scatter!(ts, xs, mc=1, shape=:circle,
      label = "h=" * string(h))

   h /= 2
end
plot!(ts, tanh.(ts), lc=2, lw = 0.5,
   label = "tanh(t)")

正確な解との誤差評価

plot(xlabel="h", xscale=:log10, yscale=:log10)

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
   scatter!([h], [e], mc=1, shape=:circle)
   h /= 2
end
xlims!(1e-2, 1)
ylims!(1e-4, 1e-1)
plot()

▶ 常微分方程式の初期値問題:修正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 Plots
plot(xlabel="t", ylabel="x")

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
scatter!(ts, xs, lc=1)
plot!(ts, tanh.(ts), lc=2)

刻みを狭くする

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

plot(xlabel="t", ylabel="x")

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
   scatter!(ts, xs, mc=1, shape=:circle,
      label = "h=" * string(h))
   h /= 2
end
plot!(ts, tanh.(ts), lc=2,
   label = "tanh(t)", lw = 0.5)
plot!()

正確な解との誤差評価

using LinearAlgebra
using Plots
plot(xlabel="h", xscale=:log10, yscale=:log10)

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
   scatter!([h], [e], mc=1, shape=:circle)
   h /= 2
end
xlims!(1e-2, 1)
ylims!(1e-5, 1e-1)
plot!()

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

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

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

using LinearAlgebra
using Plots
plot(xlabel="h", xscale=:log10, yscale=:log10)

h = 0.4
kmax = 8
hs = zeros(kmax)
e_euler = zeros(kmax)
e_meuler = zeros(kmax)

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

   # Euler
   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
   xtrue = tanh.(ts)
   e_euler[k] = norm(xs .- xtrue) / n

   # modified Euler
   xs[1] = 0 # initial condition
   for i = 1:n-1
      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
   xtrue = tanh.(ts)
   e_meuler[k] = norm(xs .- xtrue) / n
   h /= 2
end
scatter!(hs, e_euler, mc=1, shape=:circle,
   label = "Euler")
scatter!(hs, e_meuler, mc=2, shape=:circle,
   label = "modified Euler")
xlims!(1e-3, 1)
ylims!(1e-6, 1e-1)
plot!()

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

（少し難しいので，後回しにしてもよい）

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法