第9回： ▼ 総和・数値積分

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

aux = (label="", lc=:black, lw=0.5)
using LaTeXStrings
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)

▼ 級数和の公式（繰り返しで加算)

自然数の級数和の結果がいくつか知られている．これらのプロットを描いて，結果を確認しよう．

\[\sum_{k=1}^{n} k = 1 + 2 + \cdots + k + \cdots + n = \dfrac{n(n+1)}{2}\]

plot(xlabel="n", ylabel="sum i up to n")

nmax = 25
xs1 = 0:0.2:nmax
plot!(xs1, xs1 .* (xs1 .+ 1) / 2,
   label = "sum i", lc=1)

ns = 0:nmax
for n in ns
   s1 = 0.0
   for i = 1:n
      s1 += i
   end
   scatter!([n], [s1], ms=4, mc=2, shape=:circle)
end
plot!()

\[\sum_{k=1}^{n} k^2 = 1^2 + 2^2 + \cdots + k^2 + \cdots + n^2 = \dfrac{n(n+1)(2n+1)}{6}\]

plot(xlabel="n",ylabel="sum i^2 up to n")

nmax = 25
xs1 = 0:0.2:nmax
plot!(xs1, xs1 .* (xs1 .+ 1) .* (2 * xs1 .+ 1) / 6, lc=1)

ns = 0:nmax
for n in ns
   s = 0.0
   for i = 1:n
      s += i^2
   end
   scatter!([n], [s], ms=4, mc=2, shape=:circle)
end
plot!()

■ ベクトルのインデックス

Base.length — Method

参考 → ■ ベクトル

ベクトル a の寸法は，関数 length(a) で得られる．

julia> v = [11, 21, 31, 41, 51]5-element Vector{Int64}:
 11
 21
 31
 41
 51
julia> length(v)5

Indexing (section)

ベクトル a ，整数 i に対して a[i] と書くと，ベクトル a の i 番目の要素の値が得られる．要素の番号（「インデックス，index」） i は 1 から数える． end というインデックスは，ベクトルの最後の要素を指す．

julia> v[1]11
julia> v[2]21
julia> v[end] # 末尾の要素51
julia> v[end-1] # 末尾の一つ前の要素41

存在しないインデックスを指定すると，例外が発生する．

julia> v[0] # => ERROR: BoundsErrorERROR: BoundsError: attempt to access 5-element Vector{Int64} at index [0]

インデックスとして，整数 i の代わりに，範囲を指定すると，その範囲のインデックスを持つベクトルが得られる．（参考 ■ 範囲）

julia> v[2:3]2-element Vector{Int64}:
 21
 31
julia> v[1:end-1] # 最初から，末尾の一つ前の要素4-element Vector{Int64}:
 11
 21
 31
 41
julia> v[4:6] # => ERROR: BoundsErrorERROR: BoundsError: attempt to access 5-element Vector{Int64} at index [4:6]

■ ベクトルの生成

ベクトルを作る方法は，いくつかある．

これまでに，以下の方法を紹介した．

要素を列挙する方法（■ ベクトル）
範囲を用いる方法（■ 範囲）
関数 range を用いる方法．結果は範囲となる．（■ 等差数列）

■ 疑似乱数を要素とするベクトルを作る

Random.rand — Function

julia> rand(10) # => 10-elements10-element Vector{Float64}:
 0.30705996007750513
 0.5409092288457358
 0.3579872356311471
 0.9396837278185431
 0.20496424903047905
 0.2311964333481108
 0.6297967305217758
 0.16824468478622712
 0.4678969137764203
 0.683090316265765
julia> rand([1, 2, 3], 10) # [1,2,3]からランダムに10個選ぶ10-element Vector{Int64}:
 2
 2
 2
 1
 3
 3
 2
 1
 1
 1

ヒストグラムを描く．分割数 10

plot()

xs = rand(1000)
# hist!(xs, bins = 10)
plot!(xs, seriestype=:hist, bins = 10)

xlims!(-0.2, 1.2)

■ 正規乱数を要素とするベクトルを作る

Base.Random.randn — Function

平均 $0$，標準偏差 $1$ の正規分布の疑似乱数を作る

julia> randn(10) # => 10-elements10-element Vector{Float64}:
  2.953209380222828
 -0.6483941143370618
 -1.3587420013121287
 -0.3398675073494283
 -0.9362619445527262
 -0.6960672428932785
 -1.7933309864428233
  0.5007403187644065
  0.6627474226043902
 -0.08808999083627937

ヒストグラムを描く．分割数 50

plot()
xs = randn(1000)
plot!(xs, seriestype=:hist, bins = 50)
xlims!(-4, 4)

● 内包表記

julia> [x^2 for x = 0:10]11-element Vector{Int64}:
   0
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100
julia> [x^2 for x in [-3, 0, 2]]3-element Vector{Int64}:
 9
 0
 4
julia> [x^2 for x in -10:2:10 if rem(x, 3) != 2]9-element Vector{Int64}:
 100
  64
  36
  16
   4
   0
  16
  36
 100

フーリエ級数の和

▼ フーリエ級数の和（繰り返しで加算)

周期波形 $f(t+T) = f(t)$ は，以下のように，三角関数の級数和として表される．ここで，$a_0, a_1, \cdots$, $b_1, b_2, \cdots$ は実数の定数である．これを，実フーリエ級数和という．

\[\begin{aligned} f(t) & = a_0 \\ & + a_1 \cos \omega{t} + b_1 \sin \omega{t} \\ & + a_2 \cos 2\omega{t} + b_2 \sin 2\omega{t} \\ & + a_3 \cos 3\omega{t} + b_3 \sin 3\omega{t} + \cdots \end{aligned}\]

ここで $\omega$ は基本周波数である．

\[\omega=\dfrac{2\pi}{T}\]

以下の例では，すでに知られているフーリエ級数和から，元の関数が近似される様子を観察するのに留める．

▼ 方形波：フーリエ級数の有限和

方形波は，▶ 方形波を描くで紹介した．

基本周波数 $\omega=1$，数 $-1$ と $1$ とを往復する方形波を描こう．

plot()

ts = -3pi:pi/360:3pi
plot!(ts, sign.(sin.(ts)), lw=1)
yticks!([-1, 0, 1], ["-1", "0", "1"])
xticks!(
   [-3pi, -2pi, -pi, 0, pi, 2pi, 3pi],
   [L"-3\pi", L"-2\pi", L"-\pi", "0", L"\pi", L"2\pi", L"3\pi"],
)

この方形波のフーリエ級数和は，以下のように与えられる．

\[f(t) = \dfrac{4}{\pi}\left\{\sin{t}+\dfrac{\sin{3t}}{3}+\dfrac{\sin{5t}}{5}+\cdots\right\}\]

この式の $\sin t$, $\sin 3t$, $\sin 5t$ の３つを加えると，方形波に近くなることを観察する．

plot(lw=1)

ts = -3pi:pi/360:3pi
ys = sin.(ts) * 4 / pi
plot!(ts, ys, label = "n=1")
yticks!([-1, 0, 1], ["-1", "0", "1"])
xticks!(
   [-3pi, -2pi, -pi, 0, pi, 2pi, 3pi],
   [L"-3\pi", L"-2\pi", L"-\pi", "0", L"\pi", L"2\pi", L"3\pi"],
)
ys += sin.(3ts) / 3 * 4 / pi
plot!(ts, ys, label = "n=1,3")
ys += sin.(5ts) / 5 * 4 / pi
plot!(ts, ys, label = "n=1,3,5")

今度は $\sin 13t$ まで加えた結果を観察しよう．

plot(lw=1)

ts = -3pi:pi/36:3pi
n = 13
ys = zero.(ts)
for i = 1:2:n
   global ys
   ys += sin.(i * ts) / i * 4 / pi
end
plot!(ts, ys)
plot!(ts, sign.(sin.(ts)), label = "up to" * string(n))

yticks!([-1, 0, 1], ["-1", "0", "1"])
xticks!(
   [-3pi, -2pi, -pi, 0, pi, 2pi, 3pi],
   [L"-3\pi", L"-2\pi", L"-\pi", "0", L"\pi", L"2\pi", L"3\pi"],
)

Note

上のフーリエ級数和が方形波を近似すると説明したが，なめらかな三角関数の級数和をいくら加えていっても，なめらかでない方形波に一致することはない．級数和が元の関数に近づくのは「各点収束（pointwise convergence）」ではなく「一様収束（uniform convergence）」に相当する．

▼ 三角波：フーリエ級数の有限和

一定の正の傾きで増加，一定の負の傾きで減少を繰り返す周期関数を，三角波（triangular wave）という．

傾き $1$ と $-1$ で，周期 $2\pi$ の三角波を描こう．この関数は，絶対値関数 abs （参考: ▼ 絶対値関数）と関数 mod2pi （参考: ▶ 2piで割った剰余）とを組み合わせて定義できる．参考→ ■ 関数の定義（代入文形式）

triangular(t) = pi - abs.(mod2pi.(t) - pi)

plot(aspect_ratio=:equal)

ts = -3.5pi:pi/6:3.5pi
plot!(ts, triangular.(ts), lc=1)
xlims!(-pi * 2.5, pi * 2.5)
ylims!(-pi * 0.1, pi * 1.1)

上の三角波のフーリエ級数展開は，以下の通りである．

\[f(t) = \dfrac{\pi}{2} - \dfrac{4}{\pi}\left\{ \cos t + \dfrac{\cos 3t}{3^2} + \dfrac{\cos 5t}{5^2} + \cdots\right\}\]

plot(aspect_ratio=:equal)
ts = -3.5pi:pi/6:3.5pi
plot!(ts, triangular.(ts), lc=1)

ys = one.(ts) * (pi / 2)
for n = 1:2:5
   global ys
   ys -= cos.(n * ts) * (4 / pi / n^2)
end
scatter!(ts, ys, ms=4, mc=2, shape=:circle)
ylims!(-pi * 0.1, pi * 1.1)

勾配が不連続に変化する点（キンク，kink）を拡大して描画しよう．

plot(aspect_ratio=:equal)

ts = -3.5pi:pi/6:3.5pi
for nmax = 1:2:9
   ys = one.(ts) * (pi / 2)
   for n = 1:2:nmax
      ys -= cos.(n * ts) * (4 / n^2 / pi)
   end
   plot!(ts, ys, label = nmax)
end
xlims!(-pi * 0.1, pi * 2.1)
ylims!(-pi * 0.1, pi * 1.1)

◀ 練習：フーリエ級数の有限和

次の級数和で表される曲線を描け．

\[f(t) = \dfrac{4}{\pi}\left\{ \sin t - \dfrac{\sin 3t}{3^2} + \dfrac{\sin 5t}{5^2} - \cdots\right\}\]

▼ 数値積分

定積分の近似値を，級数和として求めることができる（数値積分）．

以下では，連続関数の，有限な区間に対する定積分の近似値を求めてみる．参考→ ▼ 関数が連続とは

例として，関数 $g(x)$

\[g(x) = \dfrac{1}{1+x}\]

を，$x = 0$ から $1$ の範囲で積分しよう．

関数 $g(x)$ は，この範囲で単調減少である．

plot()

xmin = 0
xmax = 1
m = 6
n = 2^m
xs = range(xmin, xmax, length = n + 1)
g(x) = 1 / (1 + x)
plot!(xs, g.(xs), lc=:black, seriestype=:steppre)
ylims!(0, 1.2)

定積分の値は，

\[\int_{0}^{1}\dfrac{1}{1+x}\;dx = \left[\log\left\vert{1+x}\right\vert\right]_{x=0}^{x=1} = \log{2}\]

である．

▼ Riemann和（繰り返しで加算)

積分すべき関数を，等間隔の短冊に区切り，短冊の面積の総和をとろう．

短冊の幅を $d$ とすると，

\[s_{1} = \sum_{i=1}^{n} g(x_{i})\cdot{d}\]

という，総和（「Riemann和」）をとることになる．

以下のプロットは，$8$ 枚の短冊に分けた様子を示す．ここで，短冊の高さは，各短冊の左端の関数の値をとった．

plot()

xmin = 0
xmax = 1
m = 3
n = 2^m
xs = range(xmin, xmax, length = n + 1)  # n個の短冊に分割する
d = (xmax - xmin) / n  # 短冊の刻み
g(x) = 1 / (1 + x)
plot!(xs, g.(xs))
ylims!(0, 1.2)

for x in xs[1:end-1]
   plot!([x, x, x + d, x + d], [0, g(x), g(x), 0], linewidth=2; aux...)
end
plot!()

では，短冊を $2^4 = 16$ 枚に分けて，短冊の面積の総和をとろう．

se = log(2)
m = 4
n = 2^m
xs = range(xmin, xmax, length = n + 1)
d = (xmax - xmin) / n
s1 = 0
for i = 1:n
   global s1
   x = xs[i]
   s1 += g(x) * d
end
#
@show s1, se, (s1 - se) / se;

(0.7090162022075267, 0.6931471805599453, 0.022894158834725318)

$16$ 分割でも，相対誤差 $2.3\%$ を達成した．

分割数を増やせば，この和は，正しい定積分の値に近づいていくであろう．

分割数を $2^{m}$ で増やして，絶対誤差を描こう．横軸の分割数は，対数で示した．

plot(xscale=:log10, xlabel="n", ylabel="absolute error")

se = log(2)
for m = 0:12
   n = 2^m
   xs = range(xmin, xmax, length = n + 1)
   d = (xmax - xmin) / n
   s1 = 0
   for i = 1:n
      x = xs[i]
      s1 += g(x) * d
   end
   scatter!([n], [abs(s1 - se)], mc= :black, shape=:circle)
end
plot!()

今度は，相対誤差を，両対数プロットで描く．

plot(xscale=:log10, yscale=:log10, xlabel="n", ylabel="relative error (absolute value)")

for m = 0:12
   n = 2^m
   xs = range(xmin, xmax, length = n + 1)
   d = (xmax - xmin) / n
   s1 = 0
   for i = 1:n
      x = xs[i]
      s1 += g(x) * d
   end
   scatter!([n], [abs(s1 - se) / se], mc= :green, shape=:circle)
end
plot!()

■ 総和関数 sum

Base.sum — Function

関数 sum(xs) は，数のコレクション $v$ を引数にとり，$v$ のすべての要素の総和を求める．

julia> sum([1, 2, 3, 4, 5])15
julia> sum(1:5)15

▼ 級数和の公式（関数 sumを用いる)

\[\sum_{k=1}^{n} k = 1 + 2 + \cdots + k + \cdots + n^2 = \dfrac{n(n+1)}{2}\]

plot(xlabel="n",ylabel="sum i up to n")

nmax = 25
xs1 = 0:0.2:nmax
plot!(xs1, xs1 .* (xs1 .+ 1) / 2, label = "sum i", lc=1)

ns = 0:nmax
for n in ns
   xs = 1:n
   s1 = sum(xs)
   scatter!([n], [s1], mc=2, shape=:circle)
end
plot!()

\[\sum_{k=1}^{n} k^2 = 1^2 + 2^2 + \cdots + k^2 + \cdots + n^2 = \dfrac{n(n+1)(2n+1)}{6}\]

plot(xlabel="n",ylabel="sum i^2 up to n")

nmax = 25
xs1 = 0:0.2:nmax
plot!(xs1,
   xs1 .* (xs1 .+ 1) .* (2 * xs1 .+ 1) / 6,
   lc=1)

ns = 0:nmax
for n in ns
   # 各要素を二乗
   xs = (1:n) .^ 2
   s = sum(xs)
   scatter!([n], [s], ms=4, mc=2, shape=:circle)
end
plot!()

▼ Riemann和（関数 sumを用いる)

Riemann和において，刻み幅 $d$ はすべての短冊に共通であるから，$d$ をくくりだして

\[s_{1} = \sum_{i=1}^{n} g(x_{i})\cdot{d} = d\cdot\sum_{i=1}^{n} g( x_{i})\]

のようにまとめることができる．すなわち，関数の値の和 $\sum_{i=1}^{n} g( x_{i})$ をとってから $d$ 倍すればよい．関数の値の和を取るのに，関数 sum を使うことができる．

下のプログラムで g.(xs[1:end-1]) は，ベクトル xs[1:end-1] の各要素に関数 g() を適用したベクトルである．

ループで和を計算した場合と，関数 sum を用いる場合との両方で，相対誤差を描く．計算結果が一致していることが見える． (参考: 結果が一致することを確かめるプロットの描画 → ●▼ 周期関数を確認する

plot(xscale=:log10, yscale=:log10,
   xlabel="n", ylabel="absolute error")

se = log(2)
for m = 0:12
   n = 2^m
   xs = range(xmin, xmax, length = n + 1)
   d = (xmax - xmin) / n

   # 和を取る
   s1 = 0
   for i = 1:n
      x = xs[i]
      s1 += g(x) * d
   end
   scatter!([n], [abs(s1 - se) / se], shape=:circle,
      ms=4, mc=1, label="summation")

   # sum を使う
   s2 = sum(g.(xs[1:end-1])) * d
   scatter!([n], [abs(s2 - se) / se], shape=:circle,
      mc=2, label="sum function")
end
plot!()

▼ 台形則（関数 sumを用いる)

今度は，短冊を台形として計算してみる．

plot()

m = 2
n = 2^m
xmin = 0
xmax = 1
xs = range(xmin, xmax, length = n + 1)
d = (xmax - xmin) / n

g(x) = 1 / (1 + x)
plot!(xs, g.(xs), lc=1)
ylims!(0, 1.2)

for i = 1:n
   x = xs[i]
   plot!([x, x, x + d, x + d], [0, g(x), g(x + d), 0]; aux...)
end
plot!()

総和をとるとき，隣り合う台形の面積をまとめることができることに注目しよう．

\[s_{t} = \sum_{i=1}^{n} \dfrac{g(x_i)+g(x_{i+1}) }{2}\cdot{d} = d\cdot\left[ \dfrac{g(x_1)}{2} + \sum_{i=2}^{n-1} g(x_i) + \dfrac{g(x_{n})}{2} \right]\]

先の Riemann和と台形則の値を両方計算してみよう．

se = log(2)

# Riemann和
s1 = 0
for x in xs[1:end-1]
   global s1
   s1 += g(x) * d
end

# 台形則
st = (g(xs[1]) + g(xs[end])) / 2
for i = 2:n
   global st
   x = xs[i]
   st += g(x)
end
st *= d
#
@show s1, st, (s1 - se) / se, (s1 - se) / se;

(0.7595238095238095, 0.6970238095238095, 0.09576123343709363, 0.09576123343709363)

関数 sum を使って簡潔に書こう．

se = log(2)
m = 4
n = 2^m
xs = range(xmin, xmax, length = n + 1)
d = (xmax - xmin) / n
s1 = sum(g.(xs[1:end-1])) * d

st = (g(xs[1]) + g(xs[end])) / 2
st += sum(g.(xs[2:end-1]))
st *= d
#
@show s1, st, (s1 - se) / se, (s1 - se) / se;

(0.7090162022075267, 0.6933912022075267, 0.022894158834725318, 0.022894158834725318)

相対誤差を描く．

plot(xscale=:log10, yscale=:log10, xlabel="n", ylabel="relative error (absolute value)")

se = log(2)
for m = 0:12
   local n = 2^m
   local xs = range(xmin, xmax, length = n + 1)
   local d = (xmax - xmin) / n

   local s1 = sum(g.(xs[1:end-1])) * d

   local st = (g(xs[1]) + g(xs[end])) / 2
   st += sum(g.(xs[2:end-1]))
   st *= d

   scatter!([n], [abs(s1 - se) / se], shape=:circle)
   scatter!([n], [abs(st - se) / se], shape=:circle)
end

◀ 練習：Riemann和・台形則

以下の定積分の近似値を，Riemann和と台形則でそれぞれ評価してみよ．（注記されていない）定積分の理論値は各自計算せよ．

\[\begin{gathered} \int^{1}_{0} 3x^2\;dx\;, \\ \int^{1}_{0} 3 \left(x+1 \right)^2\;dx\;, \\ \int^{1}_{0} \exp{x}\;dx\;, \\ \int^{2}_{0} \dfrac{1}{(1+x)^2}\;dx\;=\dfrac{2}{3}, \\ \int^{\pi}_{0} \sin{x}\;dx\;, \\ \int^{1}_{-1} \dfrac{2}{1+x^2}\;dx = \pi \end{gathered}\]

■ 繰返し内部からの脱出

Repeated Evaluation: Loops

for 文の繰り返し（「 for ブロック」）の内部で，「 break 文」を使うと，現在繰り返し中のループから直ちに抜けることができる．

julia> for i = 1:1000
          println(i)
          if i >= 5
             break
          end
       end1
2
3
4
5

乱数の値が $0.8$ を超えるまで繰り返す．

for i = 1:10
   r = rand()
   println(r)
   if r > 0.8
      break
   end
end

0.4646963625108418
0.1352433580374306
0.8563282957862134

二重ループ，内側のループからの脱出

julia> for j = 1:3
          for i = 1:5
             println("i=" * string(i) * " j=" * string(j))
             if i >= 3
                break
             end
          end
       endi=1 j=1
i=2 j=1
i=3 j=1
i=1 j=2
i=2 j=2
i=3 j=2
i=1 j=3
i=2 j=3
i=3 j=3

一つの for 文に二つの繰り返しを書いた場合，break で for 文全体から抜けてしまう．

julia> for j = 1:3, i = 1:5
          println("i=" * string(i) * " j=" * string(j))
          if i >= 3
             break
          end
       endi=1 j=1
i=2 j=1
i=3 j=1

for ブロックの内部で，continue 文を使うと，次の繰り返しに直ちに移動できる．以下で，i % 3 は rem(i,3) と同じである．参考→ ■ 残余 rem と整商 div

julia> for i = 1:10
          if i % 3 != 0
             continue
          end
          println(i)
       end3
6
9

◀● 練習：条件が成り立つまで繰り返す：数値積分

(少し難しいので，後回しにしてもよい)

分割数 $n$ を $2^{20}$ まで，$2$ の冪乗で増やしていけ，ただし，相対誤差が $10^{-4}$ 以下になったら，そこで終了せよ．

▼ Riemann和（関数 sumを用いる) ，または， ▼ 台形則（関数 sumを用いる) の，どちらを用いてもよい．

今回のまとめ

ベクトルのインデックス
要素が 0 または 1 のベクトルの生成
ベクトルの総和 sum
級数和
フーリエ級数の和
数値積分：Riemann和
数値積分：台形則
繰返し内部からの脱出

第9回： ▼ 総和・数値積分

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

▼ 級数和の公式（繰り返しで加算)

■ ベクトルのインデックス

■ ベクトルの生成

■ 要素が 0のベクトルを作る

■ 要素が 1 のベクトルを作る

■ 疑似乱数を要素とするベクトルを作る

■ 正規乱数を要素とするベクトルを作る

● 内包表記

フーリエ級数の和

▼ フーリエ級数の和（繰り返しで加算)

▼ 方形波：フーリエ級数の有限和

▼ 三角波：フーリエ級数の有限和

◀ 練習：フーリエ級数の有限和

▼ 数値積分

▼ Riemann和（繰り返しで加算)

■ 総和関数 sum

▼ 級数和の公式（関数 sumを用いる)

▼ Riemann和（関数 sumを用いる)

▼ 台形則（関数 sumを用いる)

◀ 練習：Riemann和・台形則

■ 繰返し内部からの脱出

◀● 練習：条件が成り立つまで繰り返す：数値積分

今回のまとめ