[第9回： ■ 配列要素の操作／▶常微分方程式の数値解法](@id ch09)

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

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

■ 積

Base.prod — Function

v = [2, 3, 4];
prod(v)

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

■ ノルム

LinearAlgebra.norm — Function

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

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

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

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

using LinearAlgebra
norm(v)
@show sqrt(sum(v .^ 2))

r = 0;
for i = 1:length(v)
   global r
   r += v[i]^2
end
sqrt(r)

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

v = [1, 2, 3, 4, 5, 6, 7];
abs.(v)

■ 平均値・標準偏差

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）」標準偏差を算出する．

v = [1, 2, 3, 4, 5, 6, 7];
# Statistics パッケージの読み込み
using Statistics
# 平均値
mean(v)
sum(v) / length(v)
# 偏りがない標準分散，n-1 で割る
std(v)
sqrt(sum((v .- mean(v)) .^ 2) / (length(v) - 1))
# 偏った標準分散，n で割る
std(v, corrected = false)
sqrt(sum((v .- mean(v)) .^ 2) / (length(v)))

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

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

Base.min — Function
Base.max — Function

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

■ `splatting` 演算子

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

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

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

■ ベクトル要素への代入

v = collect(1:10)
# インデックス：整数
v[4] = 0
v

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

# インデックス：範囲
v[3:2:10] .= 0
v
# `=` では例外を発生する
v[3:2:10] = 1 # => Exception

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

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

数 $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.primes` — Function](https://juliamath.github.io/Primes.jl/stable/api/#Primes.prime)
* [`Primes.isprime` — Function](https://juliamath.github.io/Primes.jl/stable/api/#Primes.isprime)

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

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

# import Pkg; Pkg.add("Primes") # コメントを外してパッケージを設置する．一度だけ行えばよい
using Primes
isprime(2)
isprime(3)
isprime(4)
isprime.([2, 3, 4])
primes(100)

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

▶ 常微分方程式の初期値問題: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

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

using PyPlot
plt.figure(); #hide

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")
plt.savefig("ch09-euler1-plot.svg"); nothing; #hide
plt.close("all") #hide

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

using PyPlot
plt.figure(); #hide

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")
plt.savefig("ch09-euler2-plot.svg"); nothing; #hide
plt.close("all") #hide

刻みを狭くする

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

using PyPlot
plt.figure(); #hide

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()
plt.savefig("ch09-euler3-plot.svg"); nothing; #hide
plt.close("all") #hide

正確な解との誤差評価

using LinearAlgebra
using PyPlot
plt.figure(); #hide

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)
plt.savefig("ch09-euler4-plot.svg"); nothing; #hide
plt.close("all") #hide

▶ 常微分方程式の初期値問題:修正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

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

using PyPlot
plt.figure(); #hide

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")
plt.savefig("ch09-meuler1-plot.svg"); nothing; #hide
plt.close("all") #hide

刻みを狭くする

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

using LinearAlgebra
using PyPlot
plt.figure(); #hide

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()
plt.savefig("ch09-meuler2-plot.svg"); nothing; #hide
plt.close("all") #hide

正確な解との誤差評価

using LinearAlgebra
using PyPlot
plt.figure(); #hide

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)
plt.savefig("ch09-meuler4-plot.svg"); nothing; #hide
plt.close("all") #hide

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

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

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

using PyPlot
plt.figure(); #hide

using LinearAlgebra
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
plt.plot(hs, e_euler, ".", label = "Euler")
plt.plot(hs, e_meuler, ".", label = "modified Euler")
plt.xlabel("h")
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.xlim(1e-3, 1)
plt.ylim(1e-6, 1e-1)
plt.savefig("ch09-meuler6-plot.svg"); nothing; #hide
plt.close("all") #hide

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

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

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法

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ch09.md

ch09.md

[第9回： ■ 配列要素の操作／▶常微分方程式の数値解法](@id ch09)

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

■ 積

■ ノルム

■ 平均値・標準偏差

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

■ `splatting` 演算子

■ ベクトル要素への代入

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

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

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

例題：Euler法による解法

刻みを狭くする

正確な解との誤差評価

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

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

刻みを狭くする

正確な解との誤差評価

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

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

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

★ 今回のまとめ

Files

ch09.md

Latest commit

History

ch09.md

File metadata and controls

[第9回： ■ 配列要素の操作／▶常微分方程式の数値解法](@id ch09)

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

■ 積

■ ノルム

■ 平均値・標準偏差

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

■ splatting 演算子

■ ベクトル要素への代入

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

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

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

例題：Euler法による解法

刻みを狭くする

正確な解との誤差評価

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

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

刻みを狭くする

正確な解との誤差評価

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

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

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

★ 今回のまとめ

■ `splatting` 演算子

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