Similar to this post I am finding that constraints tend to make FindFit poor at finding a good fit. However, unlike the previous post, the constraints are required and the success of FindFit seems to be very sensitive to initial values.
Below I describe my particular problem in detail. Any suggestions to make the parameter fitting more robust would be appreciated.
I am trying to fit the parameters A, a0, and n in the following nonlinear ODE
ExpODE = a'[t] == A * Exp[0.0583 * n * Sdot * t * Sqrt[a[t]]
to a set of experimental data. The primary varyiable a[t] is the crack length as a function of time t, and a[0]=a0. For those familiar with linear elastic fracture mechanics, this is a crack growth equation for a linearly increasing far-field stress S = Sdot * t. My experimental data consists of time to failure tf for various stress rates, Sdot. As an example, here is some fake experimental data
FakeData = {{100, 2.67886}, {200, 1.80317}, {300, 1.51754}, {400, 1.10085}, {500, 0.984215}, {600, 0.9044}, {700, 0.724567}, {800, 0.784177}, {900, 0.654773}, {1000, 0.627198}}
which looks like this when plotted
The link between my ODE and my experimental data is an equation for the crack length at failure:
a[tf] == 294.295 / (Sdot^2 * tf^2)
Given values for a0, A, and n, I can solve ExpODE using NDSolve for a[t], and solve for the time when the crack length equals the failure crack length. Here is my function for doing just that:
ExpModelFun[a0j_?NumberQ, Aj_?NumberQ, nj_?NumberQ][Sdotk_?NumberQ] :=
(Sol = NDSolve[{ExpODE /. {A -> Aj, n -> nj, Sdot -> Sdotk}, a[0] == a0j, WhenEvent[a[t] > 0.140, "StopIntegration"]}, a[t], {t, 0, 60}][[1]];
tf = t /. FindRoot[(a[t] == 294.295 / (Sdot^2 * t^2)) /. Join[Sol, {Sdot -> Sdotk}], {t, 0.01}];
ExpModelFun[a0j, Aj, nj][Sdotk] = tf)
(The WhenEvent clause is there to avoid having the crack run away on me. Once the crack length goes past 0.140 the ODE gets very stiff.) Of course, I do not know a0, A, and n a-priori, so I am using FindFit to iteratively solve for them. The following code
FindFit[FakeData, {ExpModelFun[a0, A, n][Sdot]}, {{A, 0.001}, {a0, 1*10^-9}, {n, 0.5}}, Sdot, StepMonitor :> Print[{A, a0, n}]]
results in
NDSolve::nrnum1: "The function value -0.140165-1.98689*10^-15 i is not a real number when the arguments are {1.010824975706311`*^-6,-0.000164548-1.9868921541309406`*^-15 i}"
a whole host of other errors, and a final fit of
{a0 -> -0.00013947, A -> -0.00290892, n -> 18.9563}
These parameters are not physical. I can fix this by adding the constraints a0 > 0 (initial crack length is positive) and A > 0 (cracks cannot shrink), but FindFit seems to get stuck, even if I set the AccuracyGoal to 1. The values it gets stuck at are
{a0 -> 4.58074*10^-9, A-> 0.000743962, n -> 1.37966}
The starting guesses {a0, 1*10^-9}, {A, 0.001}, {n, 0.5} are the parameters I used to create FakeData, before adding a little bit of noise, so they should be pretty good. I went ahead and changed 0.001 to 0.0005 anyways and found that
FitSol = FindFit[FakeData, {ExpModelFun[a0, A, n][Sdot], {a0 > 0, A > 0}}, {{a0, 1*10^-9}, {A, 0.0005}, {n, 0.5}}, Sdot, StepMonitor :> Print[{a0, A, n}], AccuracyGoal -> 3]
produces
{a0 -> 1.65801*10^-10, A -> 0.000756955, n -> 1.33919}
which looks like
This is a decent fit. What's frustrating is these values are quite close to the values FindFit got stuck at previously. I am guessing that FindFit is getting stuck in a local minimum. In my real application, I doubt I will have very good initial guesses. Am I just doomed to tweak initial values endlessly until I get right next to good fit?
