# Fitting data using a differential equation

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 == 650000}, y, {x, 0, 60}]

data = Table[{x, Evaluate[y[x] /. s], F[x]}, {x, 0, 60, 1}]

• @Nasser ...error corrected – thils Jan 8 '14 at 6:35
• What are the triplets in data ? Is it {x, y[x], y'[x]} ? – b.gates.you.know.what 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

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 == 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.

Update

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. You can define then a new data to compare.

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


which gives you Again, hope it helps.

• 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