comparison.jl

using PyPlot
using PyCall

PyDict(pyimport("matplotlib")["rcParams"])["font.size"] = 18

include("ode_solvers.jl")


f(x, y) = -3sin(10*x) / exp(x+y^2)

pattern = 2

if pattern == 1
    f(x, y) = x + y
    g(x) = 2e^x - x - 1
    x₀, xₑ = 0, 10
    y₀ = 1
elseif pattern == 2
    f(x, y) = y / 2x
    g(x) = sqrt(x)
    x₀, xₑ = 1, 10
    y₀ = 1
end

h = 0.01

fig = figure(figsize = (10, 10))

ax1 = fig[:add_subplot](1, 1, 1)

xs, ys_euler = euler_method(f, x₀, xₑ, y₀, h)
ax1[:plot](xs, ys_euler, label="Euler", linestyle="dotted")

xs, ys_modified_euler = modified_euler_method(f, x₀, xₑ, y₀, h)
ax1[:plot](xs, ys_modified_euler, label="Modified Euler", linestyle="dashed")

xs, ys_runge_kutta = runge_kutta_method(f, x₀, xₑ, y₀, h)
ax1[:plot](xs, ys_runge_kutta, label="Runge Kutta", linestyle="dashdot")

xs = x₀:h:xₑ
ys_ground_truth = g.(xs)
ax1[:plot](xs, ys_ground_truth, label="Ground truth", linestyle="solid")

ax1[:set_xlabel]("x")
ax1[:set_ylabel]("f(x)")

ax1[:legend]()
fig[:savefig]("comparison-partern$pattern-h$h.png")

cla()

ax2 = fig[:add_subplot](1, 1, 1)
ax2[:plot](xs, abs.(ys_ground_truth - ys_euler), label="Euler")
ax2[:plot](xs, abs.(ys_ground_truth - ys_modified_euler), label="Modified Euler")
ax2[:plot](xs, abs.(ys_ground_truth - ys_runge_kutta), label="Runge Kutta")
ax1[:set_xlabel]("x")
ax1[:set_ylabel]("error")
ax2[:legend]()
fig[:savefig]("error-partern$pattern-h$h.png")

mse(ys) = sum((ys_ground_truth - ys) .^ 2) / length(ys)

println("pattern $pattern, h = $h")
println("MSE Euler          : $(mse(ys_euler))")
println("MSE Modified Euler : $(mse(ys_modified_euler))")
println("MSE Runge Kutta    : $(mse(ys_runge_kutta))")

ode_solvers.jl

function init_xs_ys(x₀, xₑ, h)
    """
    Generate initial x values `xs` at equal interval `h` in range `[x₀, xₑ]`
    and corresponding initial y values `ys`.
    """
    xs = collect(x₀:h:xₑ)
    N = length(xs)
    ys = Vector{Float64}(N)
    xs, ys, N
end


function euler_method(f, x₀, xₑ, y₀, h)
    """
    Approximate the solution of the ordinary differential equation
    `y' = f(x, y)` with initial condition `y(x₀) = y₀` using
    Euler method, in the interval `[x₀, xₑ]` with step size `h`.
    """
    xs, ys, N = init_xs_ys(x₀, xₑ, h)
    x, y = x₀, y₀
    for i in 1:N
        ys[i] = y

        y = y + h*f(x, y)
        x += h
    end
    xs, ys
end


function modified_euler_method(f, x₀, xₑ, y₀, h)
    """
    Approximate the solution of the ordinary differential equation
    `y' = f(x, y)` with initial condition `y(x₀) = y₀` using
    modified Euler method, in the interval `[x₀, xₑ]` with step size `h`.
    """
    xs, ys, N = init_xs_ys(x₀, xₑ, h)
    x, y = x₀, y₀
    for i in 1:N
        ys[i] = y

        k₁ = h * f(x, y)
        k₂ = h * f(x + h, y + k₁)
        y = y + (k₁ + k₂) / 2
        x += h
    end
    xs, ys
end


function runge_kutta_method(f, x₀, xₑ, y₀, h)
    """
    Approximate the solution of the ordinary differential equation
    `y' = f(x, y)` with initial condition `y(x₀) = y₀` using
    Runge-Kutta method,  in the interval `[x₀, xₑ]` with step size `h`.
    """
    xs, ys, N = init_xs_ys(x₀, xₑ, h)
    x, y = x₀, y₀
    for i in 1:N
        ys[i] = y

        k₁ = f(x, y)
        k₂ = f(x + h/2, y + h*k₁/2)
        k₃ = f(x + h/2, y + h*k₂/2)
        k₄ = f(x + h, y + h*k₃)
        y = y + h * (k₁ + 2k₂ + 2k₃ + k₄) / 6
        x += h
    end
    xs, ys
end