E..2 Euler法のプログラム例 (1次元)

% cc euler1.c
% ./a.out
N=10
t=0.000000, x=1.000000
t=0.100000, x=1.100000
t=0.200000, x=1.210000
t=0.300000, x=1.331000
t=0.400000, x=1.464100
t=0.500000, x=1.610510
t=0.600000, x=1.771561
t=0.700000, x=1.948717
t=0.800000, x=2.143589
t=0.900000, x=2.357948
t=1.000000, x=2.593742
%

$N$ の値を大きくすると、 $t=1$ での値が $e^t=e^1=e=2.7182818284\cdots$ に近づくはずである。

実際、この微分方程式は簡単なので

$$x_N=\left(1+\frac{1}{N}\right)^N$$

であることが直接計算で分かり、 $\dsp\lim_{N\to\infty}x_N=e=x(1)$ が成り立つことが納得できる。

Julia言語によるプログラム例 euler1.jl

# euler1.jl --- dx/dt=x (0<t<1), x(0)=1 を Euler法で解く

#=
  入手の仕方
    curl -O http://nalab.mind.meiji.ac.jp/~mk/program/ode_for_ms/euler1.jl
  実行の仕方
  (1) julia euler1.jl
  (2) julia コマンドの中で include("euler1.jl") とする。
     (そうすると testeuler() が実行されるが)
      関数が定義されているので testeuler(1,100) のように呼び出せる。
=#

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

# デフォルトの Tmax=1.0, N=10 で実行
testeuler()

Julia言語によるプログラム例 euler2.jl

# euler2.jl --- dx/dt=x (0<t<1), x(0)=1 を Euler法で解く

#=
  入手の仕方
    curl -O http://nalab.mind.meiji.ac.jp/~mk/program/ode_for_ms/euler2.jl
  実行の仕方
  (1) julia euler2.jl
  (2) julia コマンドの中で include("euler2.jl") とする。
      (そうすると testeuler() が実行されるが)
      関数が定義されているので testeuler(1,100) のように呼び出せる。
=#

using Printf

function euler(f,t,x,dt)
  x + dt * f(t,x)
end

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 = euler(f,t,x,dt)
    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

# デフォルトの Tmax=1.0, N=10 で実行
testeuler()

ターミナルから Julia の REPL を起動してプログラムを読み込み、 実行する

% julia
julia> include("euler1.jl")
t=0.0000 x=1.0000000 exact=1.0000000 error=0.000e+00
t=0.1000 x=1.1000000 exact=1.1051709 error=5.171e-03
t=0.2000 x=1.2100000 exact=1.2214028 error=1.140e-02
t=0.3000 x=1.3310000 exact=1.3498588 error=1.886e-02
t=0.4000 x=1.4641000 exact=1.4918247 error=2.772e-02
t=0.5000 x=1.6105100 exact=1.6487213 error=3.821e-02
t=0.6000 x=1.7715610 exact=1.8221188 error=5.056e-02
t=0.7000 x=1.9487171 exact=2.0137527 error=6.504e-02
t=0.8000 x=2.1435888 exact=2.2255409 error=8.195e-02
t=0.9000 x=2.3579477 exact=2.4596031 error=1.017e-01
t=1.0000 x=2.5937425 exact=2.7182818 error=1.245e-01

julia> testeuler(1,100)
t=0.0000 x=1.0000000 exact=1.0000000 error=0.000e+00
t=0.0100 x=1.0100000 exact=1.0100502 error=5.017e-05
...
(中略)
...
t=0.9900 x=2.6780335 exact=2.6912345 error=1.320e-02
t=1.0000 x=2.7048138 exact=2.7182818 error=1.347e-02

julia>