第9回： ■ 配列要素の操作・▶常微分方程式の数値解法
■ ベクトルを引数とする関数
■ 総和関数 sum のように、ベクトルを引数とする関数がある。
■ 積
julia> v=[2,3,4];
julia> prod(v)
24
julia> r=1;
julia> for i in 1:length(v)
r *= v[i]
end
julia> r
24■ ノルム
julia> v=[1,2,3,4,5,6,7];
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)
r += v[i]^2
end
julia> sqrt(r)
11.832159566199232■ 平均値・標準偏差
julia> v=[1,2,3,4,5,6,7];
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■ 複数の数を引数とする関数
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(::AbstractArray{T1<:Real,N} where N, !Matched::Real) where T1<:Real at deprecated.jl:56
min(::AbstractArray{T1<:Real,N} where N, !Matched::AbstractArray{T2<:Real,N} where N) where {T1<:Real, T2<:Real} at deprecated.jl:56
min(::Any, !Matched::Any) at operators.jl:361
...
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
0
julia> v
10-element Array{Int64,1}:
1
2
0
0
0
6
0
8
0
10■ 素数の生成：エラトステネスの篩
エラトステネスの篩(ふるい)は、素数を算出する方法の一つである。 以下の手順による。
数$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
end2
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
end2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
97Note
Julia には、素数を高速に計算する関数を含むパッケージが用意されている。
primes(n) は、数 $n$ までの素数を計算する。
isprime(x)は、数 $x$ が素数であるかどうかを判定する。
julia> # 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法
問題
\[\begin{align*}
\dfrac{dx}{dt} & = 1-x^2, \\
x(0) & = 0, \\
& 0 \leq t \leq 1.6
\end{align*}\]
Euler 法による数値解法。
ただし、$t_1, t_2, \ldots$ は、一定間隔 $h$ とする。
\[\begin{align*}
\dfrac{dx}{dt} & =f(x,t), \\
t & = t_1, t_2, \ldots \\
\dfrac{x_{n+1}-x_{n}}{h} & = f(x_n,t_n) \\
x_{n+1} & = x_{n} + h f(x_n,t)
\end{align*}\]
まずは解いてみる
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
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
t=ts[i]
plot(t, x_now, "b.")
x_next=x_now+h*f(x_now, t)
@show t, x_next
x_now=x_next
end
plot(ts, tanh.(ts), "r")
xlabel("t")
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
xs=zeros(ts)
xs[1]=0 # initial condition
n=length(ts)
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
plot(ts, xs, ".")
plot(ts, tanh.(ts), "r")
xlabel("t")
ylabel("x")
刻みを狭くする
using PyPlot
tmin=0
tmax=1.6
h=0.4
for k in 1:4
ts=tmin:h:tmax
xs=zeros(ts)
xs[1]=0 # initial condition
n=length(ts)
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
plot(ts, xs, ".", label="h="*string(h))
h /= 2
end
xlabel("t")
ylabel("x")
plot(ts, tanh.(ts), "b", label="tanh(t)", lw=0.5)
legend()
正確な解との誤差評価
using PyPlot
tmin=0
tmax=1.6
h=0.4
for k in 1:5
ts=tmin:h:tmax
xs=zeros(ts)
xs[1]=0 # initial condition
n=length(ts)
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
plot(h,e,".")
h /= 2
end
xlabel("h")
xscale("log")
yscale("log")
xlim(1e-2,1)
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 法
\[\begin{align*}
\dfrac{dx}{dt} & =f(x,t), \\
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*}\]
解いてみる
#
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
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
xs=zeros(ts)
n=length(ts)
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
plot(ts, xs, ".")
plot(ts, tanh.(ts))
xlabel("t")
ylabel("x")
刻みを狭くする
using PyPlot
h=0.4
for k in 1:4
ts=tmin:h:tmax
xs=zeros(ts)
n=length(ts)
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
plot(ts, xs, ".", label="h="*string(h))
h /= 2
end
xlabel("t")
ylabel("x")
plot(ts, tanh.(ts), "b", label="tanh(t)", lw=0.5)
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 PyPlot
h=0.4
for k in 1:4
ts=tmin:h:tmax
xs=zeros(ts)
n=length(ts)
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
plot(h,e, ".")
h /= 2
end
xlabel("h")
xscale("log")
yscale("log")
xlim(1e-2,1)
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法