4
$\begingroup$

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

FakeData

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

fit

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?

$\endgroup$
1
  • $\begingroup$ hey, there's a couple of problems here but the two most obvious are that you are trying to solve a differential equation before what you call failure time (that's one) and the other is that even if you do, the minimizer finds multiple minima (i.e. your forward problem isn't well posed and there are infinite initial conditions that lead to the observation). I suggest you non-dimensionalise your ODE to see this. $\endgroup$
    – gpap
    Commented Jul 30, 2014 at 8:59

0