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

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

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

■ 積

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

julia> prod(v)
24

julia> r=1;

julia> for i in 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 in 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:104
  min(::Any, !Matched::Any) at operators.jl:414
  min(::Any, !Matched::Any, !Matched::Any, !Matched::Any...) at operators.jl:502
  ...

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> # インデックス：Range
       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: MethodError: no method matching setindex_shape_check(::Int64, ::Int64)
Closest candidates are:
  setindex_shape_check(!Matched::AbstractArray{#s72,1} where #s72, ::Integer) at indices.jl:218
  setindex_shape_check(!Matched::AbstractArray{#s72,1} where #s72, ::Integer, !Matched::Integer) at indices.jl:221
  setindex_shape_check(!Matched::AbstractArray{#s72,2} where #s72, ::Integer, !Matched::Integer) at indices.jl:225
  ...

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

エラトステネスの篩(ふるい)は、素数を算出する方法の一つである。 以下の手順による。

  • 数$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 in 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 in 2:nmax
  if sieve[i] > 0
    # println(i)
    sieve[i*2:i:nmax].=0
  end
end

for i in 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}:
  true
  true
 false

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{align*} \dfrac{x_{n+1}-x_{n}}{h} & = f(x_n,t_n) \\ x_{n+1} & = x_{n} + h f(x_n,t_n) \end{align*}

例題：Euler法による解法

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

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

刻み $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 in 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 in 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 in 1:n-1
  global x_now
  t=ts[i]
  x_now=xs[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 in 1:4
  global h
  ts=tmin:h:tmax
  n=length(ts)
  xs=zeros(n)
  xs[1]=0 #  initial condition

  for i in 1:n-1
    t=ts[i]
    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 in 1:5
  global h
  ts=tmin:h:tmax
  n=length(ts)
  xs=zeros(n)
  xs[1]=0 #  initial condition

  for i in 1:n-1
    t=ts[i]
    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{align*} 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{align*}

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

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

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

刻み $h = 0.4$ とする。

#
tmin=0
tmax=1.6
h=0.4
ts=tmin:h:tmax
x_now=0 # initial condition

n=length(ts)
for  i in 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 in 1:n-1
  global x_now, xs
  t=ts[i]
  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
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 in 1:4
  global h
  ts=tmin:h:tmax
  n=length(ts)
  xs=zeros(n)
  xs[1]=0 # initial condition
  for  i in 1:n-1
    t=ts[i]
    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 in 1:4
  global h
  ts=tmin:h:tmax
  n=length(ts)
  xs=zeros(n)
  xs[1]=0 # initial condition
  for  i in 1:n-1
    t=ts[i]
    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{align*} \dfrac{dx}{dt} & = x^2, \\ x(0) & = \dfrac{1}{2}, \\ 0 & \le t < 2 \end{align*}

解析解は、

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

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

★ 今回のまとめ

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