nonlinear-fit.jl

nonlinear-fit.jl

# Julia Version 0.4.0-dev+5181
# Commit 98c144b (2015-06-03 14:28 UTC)
#
#  - PyPlot                        1.5.3+             master
#  - PyCall                        0.8.1+             master

using Calculus, PyPlot, Roots

# function being approximated
# still needing to find a way to remove this redundancy..
cv(T::Float64, θ::Float64) = 3*k_B * (θ/T)^2 * exp(θ/T) / (exp(θ/T)-1)^2

# exchange a symbol on a expression to a given value
# general implementation using dictionaries is also possible
# credit to @ptb: http://stackoverflow.com/a/29780738/3667015
exreplace!(ex, s, v) = ex
exreplace!(ex::Symbol, s, v) = s == ex ? v : ex
function exreplace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = exreplace!(ex.args[i], s, v)
    end
    ex
end
exreplace(ex, s, v) = exreplace!(copy(ex), s, v)

# data supplied from the experimentation
k_B   = 1.9872041                                             # cal/mol*K
temps = [222.4, 283.7, 331.3, 358.5, 413.0, 520.0, 1079.7]    # kelvin
cps   = [0.762, 1.354, 1.838, 2.118, 2.661, 3.631, 5.387]     # cal/mol*K

trange = 100:0.1:3100                          # uniform distribution prior

# s(θ) = sum (y_k - c_v(T_k; theta))^2
# and ds is the derivative of s w/ respect to θ
ds_ex = :(0.)
for i in 1:7
    T = temps[i]; y = cps[i];
    ds_ex = :($ds_ex + ($y - 3*$k_B * (θ/$T)^2 * exp(θ/$T) / (exp(θ/$T)-1)^2)^2)
end
ds_ex = differentiate(ds_ex, :θ)
# probably not the best implementation, reeeeally slow here..
dS(value::Float64) = eval(exreplace(ds_ex, :θ, value))

# discover saddle points on [100, 3100]
@time saddles = fzeros(dS, 100, 3100)          # 34.308 seconds      (14645 k allocations: 554 MB, 0.32% gc time)

# mmq residue function: s(θ) = sum (y_k - c_v(T_k; theta))^2
s_ex = :(0.)
for i in 1:7
    T = temps[i]; y = cps[i];
    s_ex = :($s_ex + ($y - 3*$k_B * (θ/$T)^2 * exp(θ/$T) / (exp(θ/$T)-1)^2)^2)
end
s(value::Float64) = eval(exreplace(s_ex, :θ, value))

# check which saddle point is in fact the minimun
opt = Inf
for i in saddles
    if s(i) < s(opt) || isnan(s(opt))
        opt = i
    end
end

@time real_cps = Float64[cv(i, opt) for i in trange]  # 55.348 milliseconds (538 k allocations: 10992 KB, 15.34% gc time)

# Nonlinear Approximation Plot
axhline(y=3*k_B, linewidth=1, color="green", linestyle="--", label=L"$3k_B$ limit")
plot(temps, cps,  color="red", marker="o", linestyle="None", label=L"Measured $c_p$")
plot(collect(trange), real_cps,  color="blue", linestyle="-", label=L"Theoretic $c_p$")
title("Heat Capacity Approximation")
xlim(xmax=3100)
ylim(ymax=7)
xlabel(L"Temperatures $T$")
ylabel(L"Heat Capacities $c_p$")
legend(loc="lower right")
savefig("approx_theta_$opt.png")

@time dss      = Float64[dS(i)      for i in trange]  # 276.458 seconds      (108 M allocations: 4030 MB, 0.33% gc time)

# Verification of the dS values
axhline(y=0, linewidth=1, color="green", linestyle="--", label=L"$y=0$")
axvline(x=opt, ymin=-2, ymax=2, color="blue", linestyle="--", label="Optimal Value")
plot(collect(trange), dss, color="red", linestyle="-", label="Derivative of the MMQ")
title("MMQ Plot")
xlim(100, 3100)
xlabel(L"$\theta$ Values")
ylabel(L"$\frac{dS}{d\theta}$ Values")
legend(loc="lower right")
savefig("verifymmq_theta_$opt.png")
close()