B..1 よくある課題: Euler法、Runge-Kutta法で $x'=x$, $x(0)=1$ を解け

$$\frac{\D x}{\D t}=x(t)$$   ($t\in[0,T]$)$$,\quad x(0)=1$$

の解は

$$x(t)=e^t.$$

微分方程式の初期値問題

$$\frac{\D x}{\D t}=f(t,x) t\in[0,T] ,$$

$$x(0)=x_0$$

に対する前進 Euler 法とは

$$\Delta t:=\frac{T}{N},\quad t_j=j\Delta t$$   ( $j=0,1,\dots,N$)$$,$$

$$x_{j+1}=x_j+\Delta t f(t_j,x_j)$$   ( $j=0,1,\dots,N-1$)$$.$$

testeuler.jl

# testeuler.jl --- dx/dt=x (0<t<1), x(0)=1 を Euler法で解く
# using Printf が必要

function testeuler(Tmax=1.0,N=10)
  # 初期値
  t=0.0
  x=1.0
  @printf "t=%.4f x=%.7f exact=%.7f error=%.3e\n" t x exp(t) abs(exp(t)-x)
  dt=Tmax/N
  for j=1:N
    x += dt * f(t,x)
    t = j*dt
    @printf "t=%.4f x=%.7f exact=%.7f error=%.3e\n" t x exp(t) abs(exp(t)-x)
  end
end

function f(t,x)
  x
end

実行例

julia> using Printf
julia> include("testeuler.jl")        
julia> testeuler(1,10)
julia> testeuler(1,100)

実は $x_1=1+\frac{1}{N}$, $x_N=\left(1+\frac{1}{N}\right)^N$ となっている。 $x_N$ が $x(t_N)=x(1)=e^t=e=2.7182818284\cdots$ に近いことが確かめられる。

Euler 法の次数は $1$ であり、 誤差 $x(1)-x_N=O\left(\frac{1}{N}\right)$ であることも知られていて、 $N$ を変えることで確かめられる ($N$ を $10$ 倍にすると誤差が $\frac{1} {10}$ 倍になる)。

(4段4次)Runge-Kutta 法は

$$k_1:=\Delta t f(t_j,x_j),$$

$$k_2:=\Delta t f(t_j+\Delta t/2,x_j+k_1/2),$$

$$k_3:=\Delta t f(t_j+\Delta t/2,x_j+k_2/2),$$

$$k_4:=\Delta t f(t_j+\Delta t,x_j+k_3),$$

$$x_{j+1}:=x_j+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right).$$

testrungekutta.jl

# testrungekutta.jl --- dx/dt=x (0<t<1), x(0)=1 を Runge-Kutta法で解く
# using Printf が必要

function testrungekutta(Tmax=1.0,N=10)
  # 初期値
  t=0.0
  x=1.0
  @printf "t=%.4f x=%.7f exact=%.7f error=%.3e\n" t x exp(t) abs(exp(t)-x)
  # Runge-Kutta法
  dt=Tmax/N
  for j=1:N
    k1=dt*f(t,x)
    k2=dt*f(t+dt/2, x+k1/2)
    k3=dt*f(t+dt/2, x+k2/2)
    k4=dt*f(t+dt, x+k3)
    x += (k1+2*k2+2*k3+k4) / 6
    t = j * dt
    @printf "t=%.4f x=%.7f exact=%.7f error=%.3e\n" t x exp(t) abs(exp(t)-x)
  end
end

function f(t,x)
  x
end

実行例

julia> using Printf
julia> include("testrungekutta.jl")        
julia> testrungekutta(1,10)
julia> testrungekutta(1,100)

実は $x_1=1+\frac{1}{N}+\frac{1}{2}\frac{1}{N^2}+\frac{1}{6}\frac{1}{N^3}+ \frac{1}{24}\frac{1}{N^4}$ となっている。 誤差 $x(1)-x_N=O\left(\frac{1}{N^4}\right)$ であることが知られている。 $N$ を$10$倍にすると誤差はほぼ $\dfrac{1}{10000}$ 倍になる、 ということである。