Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have data generated from two functions given below. How can I go about finding a fit of the form:

y'[x] == -(a + b/(c + d  y[x] + e F[x] ))

where a,b,c,d,e are to be determined?

F[x_] := If[x < 37.5, 2.80 + 0.0036 x, 117.86 - 6.314 x + 0.0865 x^2]

s = NDSolve[{y'[x] == -(0.0595 + 3716/(11780 + y[x])) y[x], 
   y[0] == 650000}, y, {x, 0, 60}]

data = Table[{x, Evaluate[y[x] /. s], F[x]}, {x, 0, 60, 1}]
share|improve this question
@Nasser ...error corrected – thils Jan 8 '14 at 6:35
What are the triplets in data ? Is it {x, y[x], y'[x]} ? – b.gatessucks Jan 8 '14 at 8:12
@b.gatessucks It is {x,y[x],F[x]} – thils Jan 8 '14 at 8:26
@thils, I believe you have to approximate y'[x] based on the data using finite difference coefficients. Have you considered this? – carlosayam Jan 9 '14 at 9:44
@thils, I noticed that the solution from NDSolve was differentiable already; hence my answer :) – carlosayam Jan 13 '14 at 11:32
up vote 1 down vote accepted

This seems to be a task for FindFit. First, get the solution from NDSolve as a function; let me call it solY.

F[x_] := If[x < 37.5, 2.80 + 0.0036 x, 117.86 - 6.314 x + 0.0865 x^2];

solution = 
  NDSolve[{y'[x] == -(0.0595 + 3716/(11780 + y[x])) y[x], 
    y[0] == 650000}, y, {x, 0, 60}];

solY = y /. First[First[solution]];

Now, generate the data in a format for FindFit. In this case, as you want y'[x] == -(a + b/(c + d y[x] + e F[x] )), my suggestion is to have two independent values, y[x] and F[x] giving the dependent value y'[x]. x is not required. Note that you can use the derivative solY'!

This data is generated by the following

data = Table[{solY[x], F[x], solY'[x]}, {x, 0, 60, 1}];

Now, feed that into FindFit, like

FindFit[data, -(a + b/(c + d  v1 + e v2)), {a, b, c, d, e}, {v1, v2}]

and you get

{a -> 11425.8, b -> -2.63417*10^9, c -> 1.01206*10^7, d -> -98.9974, e -> -213696.}

The derivative solY' comes from an approximation (numerical solution to NDSolve), so I would double check results by other means. Also, your parameters are non-linear and "In the nonlinear case, it finds in general only a locally optimal fit." (see documentation)

Hope it helps.


As mentioned by OP, FindFit does not find a good answer. F is itself problematic and certainly one would like to explore the data. Here some ideas on how to do this, but consider that the parameter space is 5D - so visualization is already hard.

You can use the solution to define a new function. This is a way to do it. Say you capture the solution into resp, like resp = FindFit[...];. Then you can do

g[v1_, v2_] := Evaluate[-(a + b/(c + d v1 + e v2)) /. resp];

Using ?q actually gives g with values from the solution rules.

Mathematica graphics

You can define then a new data to compare.

data2 = Map[({#1, #2, g[#1, #2]} &) @@ # &, data];
    ColorFunction -> (Directive[Red, Opacity[0.7]] &)],
    ColorFunction -> (Directive[Blue, Opacity[0.7]] &)],
PlotRange -> All]

which gives you

Mathematica graphics

Again, hope it helps.

share|improve this answer
Appear to work fine, except that I got "Failed to converge to the requested accuracy or precision within 100 iterations"...with slightly different set of a,b,c,d,e. Possibly I need to change the model slightly? – thils Jan 14 '14 at 5:20
@thils, I am not sure about your problem space. I added few ideas you can use. Good luck! – carlosayam Jan 14 '14 at 8:26

Your Answer


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

Not the answer you're looking for? Browse other questions tagged or ask your own question.