B..2 分けて書いてみる

計算部分を独立した関数とするように分けて書くことも可能である。

Euler法の場合は例えば次のようになる。

testeuler2.jl

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

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

function testeuler2(Tmax=1.0,N=10)
  # 初期値
  t=0.0
  x=1.0
  # Euler法
  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

実行例

julia> include("testeuler2.jl")
julia> testeuler2(1.0,10)

Runge-Kutta法の場合。

testrungekutta2.jl

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

# Runge-Kutta 法の1ステップ
function rungekutta(f,t,x,dt)
  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
end

function testrungekutta2(Tmax=1.0,N=10)
  # 初期値
  t=0
  x=1
  @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
    x = rungekutta(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

実行例

julia> include("testrungekutta2.jl")
julia> testrungekutta2(10)