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.

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 at 6:35
    
What are the triplets in data ? Is it {x, y[x], y'[x]} ? –  b.gatessucks Jan 8 at 8:12
    
@b.gatessucks It is {x,y[x],F[x]} –  thils Jan 8 at 8:26
1  
@thils, I believe you have to approximate y'[x] based on the data using finite difference coefficients. Have you considered this? –  caya Jan 9 at 9:44
1  
@thils, I noticed that the solution from NDSolve was differentiable already; hence my answer :) –  caya Jan 13 at 11:32

1 Answer 1

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.

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.

Mathematica graphics

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

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 at 5:20
1  
@thils, I am not sure about your problem space. I added few ideas you can use. Good luck! –  caya Jan 14 at 8:26

Your Answer

 
discard

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.