Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.