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 am using MMA V9.0

What do you suggest for this numerical example?

I am not getting good results trying to use my own procedure to find a "best" fit using the 1-norm. The 1-norm example for FindFit that uses a fitting function to fit the first 20 primes returns {1.1478,1.00474,1.05178} for the three parameters {a,b,c}.

My procedure below, formally equivalent to the fitting function in the example, returns {1.06472, -0.0518796, 1.4054}, although the fitting curves do not differ significantly.

The FindFit example finds a minimum of 16.2208 while my procedure finds 17.1524, probably not significant.

dat = Table[Prime[u], {u, 20}]

Res[pars_] := Module[{sol, y, z, mdl, res},
  sol = Simplify[
     DSolve[{y''[x] == a*c*(2*b + c*x)/(b + c*x)^2, y[0] == 0, 
       y'[0] == a*Log[b]}, y[x], x]][[1]];
  y = y[x] /. sol;
  z = y /. {a -> pars[[1]], b -> pars[[2]], c -> pars[[3]]};
  mdl = Table[z /. x -> i, {i, 1, 20}];
  res = dat - mdl;
  Norm[res, 1]

sol = NMinimize[Res[{a, b, c}], {a, b, c}]

However, I am doing this so that I can use a nonlinear ODE and do everything numerically. The function below ResI, the numerical equivalent of Res, does what it is supposed to do.

ResI[pars_] := Module[{sol, YFcn, mdl, res},
  sol = NDSolve[{y''[x] == 
        pars[[3]]*(2*pars[[2]] + 
           pars[[3]]*x)/(pars[[2]] + pars[[3]]*x)^2, y[0] == 0, 
      y'[0] == pars[[1]]*Log[pars[[2]]]}, y, {x, .5, 20.5}][[1]];
  YFcn[x_] := y[x] /. sol;
  mdl = Table[YFcn[i], {i, 1, 20}];
  res = dat - mdl;
  Norm[res, 1]

The problem I am having is that NMinimize is not picking up numerical values for the integration. The following does not work.

solI = NMinimize[ResI[{a, b, c}], {a, b, c}]
share|improve this question
Change the ResI definition to ClearAll[ResI];ResI[pars_?(VectorQ[#, NumericQ] &)] := etc and the NMinimize line to something like solI = NMinimize[ResI[{a, b, c}], {{a, 0.5, 1.4}, {b, 0.5, 1.4}, {c, 0.5, 1.4}}]. As to the reason for the first step see this answer. – Sjoerd C. de Vries May 13 '14 at 21:25
Thanks so much. Now NMinimize returns {1.25006,1.14189,0.810974} and a 1-norm minimum of 16.148. Also helpful, your answer gives me a lot to think about. – Yuri May 14 '14 at 18:38

Your Answer


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

Browse other questions tagged or ask your own question.