# Problems with LsqFit on an ODE-based model

**URL:** <https://discourse.julialang.org/t/problems-with-lsqfit-on-an-ode-based-model/91642>\
**Category:** Optimization (Mathematical)\
**Tags:** lsqfit\
**Created:** [December 14, 2022, 2:15pm UTC](https://discourse.julialang.org/t/problems-with-lsqfit-on-an-ode-based-model/91642 "2022-12-14T14:15:40Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [December 14, 2022, 2:15pm UTC](https://discourse.julialang.org/t/problems-with-lsqfit-on-an-ode-based-model/91642/1 "2022-12-14T14:15:40Z")

</div>

I have a 1D experimental dataset and an ODE-based model. The model is able to describe experimental data reasonably well by manual fitting of the parameter:  
 ![plot_95small](https://global.discourse-cdn.com/julialang/original/3X/3/8/388a9bb5588827f3289f5bbdff083b203cf2c4a5.png)

The next step now is using `LsqFit`, at first trying to fit just one parameter:

```julia
julia> fit1 = curve_fit(model, tdata, ydata, [p0]); 
julia> fit1.converged
true
julia> plot(fit1.resid)

```

![plot_15-small-residual](https://global.discourse-cdn.com/julialang/original/3X/4/8/48d0af840b8743bbeabb2e4f90e93bf41e31215e.png)

However it apparently doesn’t even try to fit:

```julia
julia> fit1.param[1] == p0
true

```

Now let’s change the initial value

```julia
julia> fit2 = curve_fit(model, tdata, ydata, [p0*1.5]);
julia> fit2.converged
true
julia> fit2.param[1] == p0*1.5
true # !!
julia> plot(fit2.resid)

```

![plot_16-big-residual](https://global.discourse-cdn.com/julialang/original/3X/c/3/c373e70bef03cc586f3a3984d64c29f6f11b1b29.png)

The Jacobians are in both cases zero:

```julia
julia> all(fit1.jacobian .== 0)
true
julia> all(fit2.jacobian .== 0)
true

```

Now I’m completely at loss ☹

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 14, 2022, 2:33pm UTC](https://discourse.julialang.org/t/problems-with-lsqfit-on-an-ode-based-model/91642/2 "2022-12-14T14:33:02Z")

</div>

Don’t use LsqFit.jl. We actively recommend against it because it’s generally not as stable as using Optimization.jl-based cost functions with a full optimizer.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [December 14, 2022, 5:48pm UTC](https://discourse.julialang.org/t/problems-with-lsqfit-on-an-ode-based-model/91642/3 "2022-12-14T17:48:06Z")

</div>

Tomorrow is a project discussion with my boss, and it would be much better if we have this way of getting our material properties in addition/comparison to those used till now. Thus I may follow the advice and switch to another package later on, but today just trying to get results, quick and dirty.

I don’t need high precision, nor even precision estimation. Providing my own numeric differentiation for the Jacobian worked: Fit now converges to the same value for any reasonable initial parameters. Still would be nice to understand, why it didn’t work without it.

```julia
function j_m(t,p)
    J = Array{Float64,2}(undef, length(t),length(p))
    δrel = 1e-4
    p0 = p[1]
    δpm = p0 * δrel/2
    J[:,1] = (model(t, [p0+δpm]) .- model(t, [p0-δpm])) ./ (2δpm)   
    return J
end

```