# Using FindFit to fit nonlinear ODE parameters, subject to a constraint (an inverse problem)

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?

-