Hi everyone,
It’s a long shot as the problem is quite niche, but I would still appreciate any help/comments from the community.
Problem Description
To put it briefly, I am interested in fitting a nonlinear model S(\lambda; C^v; \boldsymbol{\theta}) to data pairs (S_i,\, \lambda_i). The other variable C^v is the solution of a nonlinear ODE
where \boldsymbol{\theta} are the parameters I am eventually interested in. It is basically a constitutive model from nonlinear solid mechanics. The most natural way is to integrate the ODE numerically, get C^v and then solve a nonlinear least-squares problem. I have hard-coded my own time-integrator (a 5th order RK) and it seems to work fine (I know \boldsymbol{\theta} beforehand and just ran the forward model to see if it was working fine). Using DifferentialEquations.jl is next on my list. For the fitting, I am trying to use LsqFit.jl and specifically it’s curve_fit function along with the RK5 integrator.
Data
1.000000000000000000e+00,0.000000000000000000e+00
1.080078125000000000e+00,1.007812500000000000e+01
1.157226562500000000e+00,1.660937500000000000e+01
1.255859375000000000e+00,2.221875000000000000e+01
1.336914062500000000e+00,2.584375000000000000e+01
1.452148437500000000e+00,2.982812500000000000e+01
1.618164062500000000e+00,3.306250000000000000e+01
1.771484375000000000e+00,3.612500000000000000e+01
1.945312500000000000e+00,3.881250000000000000e+01
2.111328125000000000e+00,4.153125000000000000e+01
2.302734375000000000e+00,4.421875000000000000e+01
2.515625000000000000e+00,4.743750000000000000e+01
2.716796875000000000e+00,5.050000000000000000e+01
2.996093750000000000e+00,5.446875000000000000e+01
2.923828125000000000e+00,5.156250000000000000e+01
2.800781250000000000e+00,4.721875000000000000e+01
2.664062500000000000e+00,4.343750000000000000e+01
2.535156250000000000e+00,3.962500000000000000e+01
2.400390625000000000e+00,3.584375000000000000e+01
2.259765625000000000e+00,3.168750000000000000e+01
2.156250000000000000e+00,2.879687500000000000e+01
2.054687500000000000e+00,2.553125000000000000e+01
1.956054687500000000e+00,2.118750000000000000e+01
1.845703125000000000e+00,1.757812500000000000e+01
1.743164062500000000e+00,1.304687500000000000e+01
1.584960937500000000e+00,6.894531250000000000e+00
1.469726562500000000e+00,7.329101562500000000e-01
Code (S11 is the function)
using LsqFit;
using Plots;
using DelimitedFiles;
using Optim;
function S11(λ, Cv11, params)
mu1, mu2, alph1, alph2, m1, m2, a1, a2, K1, K2, bta1, bta2, eta0, etaInf = params
l2 = λ^(-0.5)
I1 = λ^2. + 2/λ
Cv22 = Cv11^(-0.5)
Ie1 = λ^2/Cv11 + 2*l2^2/Cv22
p = ((mu1 * (I1/3.)^(alph1 - 1.) + mu2 * (I1/3.)^(alph2 - 1.)) * l2 + (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.)) * l2/Cv22) * l2
s1 = ((mu1 * (I1/3.)^(alph1 - 1.) + mu2 * (I1/3.)^(alph2 - 1.)) * λ + (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.)) * λ/Cv11) - p/λ
return s1;
end
function Cvdot11(λ, Cv11, params)
mu1, mu2, alph1, alph2, m1, m2, a1, a2, K1, K2, bta1, bta2, eta0, etaInf = params
@assert λ > 0 "Stretches can't be negative"
l2 = λ^(-0.5)
# @assert Cv11 > 0 "Cv can't be negative" # this is where it fails
Cv22 = Cv11^(-0.5)
Ie1 = λ ^ 2. /Cv11 + 2.0 * l2^2. /Cv22
A2 = (m1 * (Ie1/3.)^(a1 - 1.) + m2 * (Ie1/3.)^(a2 - 1.))
Iv1 = Cv11 + 2. * Cv22
Ie1 = λ^2/Cv11 + 2*l2^2/Cv22
Cvinv_C_sq = λ^4/Cv11^2 + 2. * l2 ^4/Cv22^2
Ie2 = 1/2. * (Ie1 ^ 2. - Cvinv_C_sq)
J2NEq = (Ie1^2/3. - Ie2) * A2^2.
etaK = etaInf + (eta0 - etaInf + K1 * (Iv1 ^ (bta1) - 3^(bta1)))/(1. + (K2 * J2NEq)^(bta2))
cvd11 = A2/etaK * (λ^2. - 1/3. * Ie1 * Cv11)
return cvd11;
end
function kterms(dt, Cv_n, λ, λ_n, λ_half, λ_quarter, λ_three_quarter, params)
G1 = Cvdot11(λ_n, Cv_n, params)
# @assert G1 > 0.
Cv_G2 = Cv_n + G1 * dt / 2.0
G2 = Cvdot11(
λ_half,
Cv_G2, params)
# @assert G2 > 1.
Cv_G3 = Cv_n + (3.0 * G1 + G2) * dt / 16.0
G3 = Cvdot11(
λ_quarter,
Cv_G3,
params)
# @assert G3 > 1.
Cv_G4 = Cv_n + G3 * dt / 2.0
G4 = Cvdot11(
λ_half,
Cv_G4,
params)
# @assert G4 > 1.
Cv_G5 = Cv_n + 3.0 * dt / 16.0 * (-G2 + 2.0 * G3 + 3.0 * G4)
G5 = Cvdot11(
λ_three_quarter,
Cv_G5,
params)
# @assert G5 > 1.
Cv_G6 = Cv_n + dt / 7.0 * (G1 + 4.0 * G2 + 6.0 * G3 - 12.0 * G4 + 8.0 * G5)
G6 = Cvdot11(
λ,
Cv_G6,
params)
# @assert G6 > 1.
Cv = Cv_n + dt / 90.0 * (7.0 * G1 + 32.0 * G3 + 12.0 * G4 + 32.0 * G5 + 7.0 * G6)
return Cv;
end
function correct_params()
params = [
(13.54 * 10^3),
(1.08 * 10^3),
(1.),
(-2.474),
(5.42 * 10^3),
(20.78 * 10^3),
(-10.),
(1.948),
(3507 * 10^3),
(10^(-6)),
(1.852),
(0.26),
(7014 * 10^3),
(0.1 * 10^3)];
return params;
end
function fit(λvals, params)
λdot = 0.01;
S11v = zeros(size(λvals, 1));
Cvvals = ones(size(λvals, 1));
for i ∈ range(1, size(λvals, 1) - 1)
Δt = abs(λvals[i+1] - λvals[i])/λdot;
λn = λvals[i];
λnp1 = λvals[i+1];
l1half = λn + 0.5 * (λnp1 - λn);
l1quarter = λn + 0.25 * (λnp1 - λn)
l1three_quarter = λn + 0.75 * (λnp1 - λn)
Cvvals[i+1] = kterms(Δt, Cvvals[i], λnp1, λn, l1half, l1quarter, l1three_quarter, params)
S11v[i+1] = S11(λnp1, Cvvals[i+1], params)
end
return S11v;
end
data = readdlm("VHBdataJulia.txt", ',');
function plot_correct_plot(data)
λdot = 0.01;
timevals = [0.0];
for i ∈ 1:length(data[:, 1])-1
Δt = abs(data[i+1, 1] - data[i, 1])/λdot;
push!(timevals, timevals[i] + Δt);
end
length(timevals)
size(data[:, 1])
plot(timevals, data[:, 1]);
savefig("time.png");
size(data[:, 1])
# trial figure
params = correct_params();
Δtf = 0.01;
λdotf = 0.01;
λfinal = 3.0;
Tf = (λfinal - 1.0)/λdotf;
nsteps = Tf/Δtf + 1;
λvals = LinRange(1.0, λfinal, Int(nsteps));
λvals = vcat(λvals, LinRange(λfinal-λdotf * Δtf, 1.0, Int(nsteps)));
λvals
plot(λvals)
savefig("λvals.png");
plot(λvals, fit(λvals, params)./1E3, label="Correct Parameters");
scatter!(data[:, 1], data[:, 2], label="Data", legend=:bottomright);
savefig("correct_params.png");
end
plot_correct_plot(data);
p0 = [13E3, 1E3, 0.9, -3., 5E3, 20E3, -12., -10., 5E6, 1.5E-6, 2., 1., 1E7, 5E1];
ib = [1E2, 1E2, 0.5, -3.5, 1E2, 1E2, -15., -15., 1E6, 1E-7, 0.5, 0.1, 1E6, 1E1];
ub = [2E4, 2E4, 1.5, 2.5, 1E5, 1E5, 3., 3., 1E8, 3E-6, 2.5, 1.0, 1E7, 2E2];
myfit = curve_fit(fit, data[:, 1], data[:, 2] .* 1E3, p0, lower=ib, upper=ub);
Result
If I run plot_correct_data above, nothing goes wrong (which led me to believe that the time-integrator was OK). Because I am calling the same function (which does the time-integration) but with the right set of parameters.
Error
But If I run curve_fit on the data, it seems to throw an error specifically at the assertion statement I had (I was getting a DomainError before and I traced it back to the second @assert statement inside the Cvdot11 function). Here is the earlier stacktrace
ERROR: LoadError: DomainError with -0.30365595683307434:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(x::Float64)
@ Base.Math .\math.jl:37
[2] ^
@ .\math.jl:909 [inlined]
[3] Cvdot11(λ::Float64, Cv11::Float64, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:22
[4] kterms(dt::Float64, Cv_n::Float64, λ::Float64, λ_n::Float64, λ_half::Float64, λ_quarter::Float64, λ_three_quarter::Float64, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:70
[5] fit(λvals::Vector{Float64}, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:110
Now, if I try to do something like (Cv + 0im)^(-0.5) I get an InexactError, namely
ERROR: LoadError: InexactError: Float64(0.5141891571782113 - 1.6389735916551083e-6im)
Stacktrace:
[1] Real
@ .\complex.jl:44 [inlined]
[2] convert
@ .\number.jl:7 [inlined]
[3] setindex!
@ .\array.jl:903 [inlined]
[4] fit(λvals::Vector{Float64}, params::Vector{Float64})
@ Main d:\curveFit\curveFit.jl:112
which seems really weird since line 112 is
Cvvals[i+1] = kterms(Δt, Cvvals[i], λnp1, λn, l1half, l1quarter, l1three_quarter, params).
It would be great if you could point me to what could be going wrong and where. Thanks in advance!
