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 2 differential equations with 2 variables, x and y,which are a function of t and I have the parameters k1, k2 y k3.

                       dx/dt=-k1 x2+ k2 x y
                      dy/dt=k1 x2-k2 x y- k3 y

I have to adjust the equations to the following experimental data

                       xo=70.26, x(t=720)=45.78
                       xo=71.04, x(t=720)=46.32
                       xo=37.23, x(t=720)=24.67
                       xo=37.91, x(t=720)=28.78

I tried FindFit and NMinimize. The problem with FindFit is that I have multiple initial conditions. Then I used NMinimize, I tried to create a function error only with the first data (eventually, I will use the rest of the data) but NMinimize gives me the following error This is the funtion.

Remove["Global`*"] 
f[k1_?NumericQ,k2_?NumericQ,k3_?NumericQ]:=
     Module[{x,y,out},out=Abs[45.78-x[720] /.
           NDSolve[{x'[t]==-k1 x[t]^2+k2 x[t]y[t],y'[t]==k1 x[t]^2-k2 x[t]y[t]-k3 y[t],
                    x[0]==70.26,y[0]==0},{x,y},{t,0,800}]]]
res=NMinimize[{f},{k1,k2,k3}]

And this is the error

NMinimize::nnum: The function value f is not a number at {k1,k2,k3} = {8.17269,8.09533,4.9417}. >>

If anybody can help

share|improve this question
    
ParametricNDSolve[] is what seems to be used now, if you have version 9. –  J. M. May 11 '13 at 2:08
    
Might the problem be caused by providing NMinimize with the list {f} as argument instead of just f? –  Sjoerd C. de Vries May 11 '13 at 6:55
    
First, the output form of NDSolve is nested list (Such as: {{x->InterpolatingFunction[{{0.,30.}},<>]}}). So you need to add such as First@, to the NDSolve to get value not {value}. Second, you need to write the function explicitly(i.e.f[k1,k2,k3] not f,) in the NMinimize. Although after these two steps, you'll still get a bunch of warnings. –  luyuwuli May 11 '13 at 10:07
    
I tried by using ParametricNDSolve and FindRoot on simpler situations; It works. But I can't find appropriate initial values of k1,k2,k3 (and y0), so Mathematica returns a bunch of warning too. I think a good constrains of k1,k2,k3 and y0 is needed. (BTW, the original equation is wrong with x^2 not x2) –  luyuwuli May 11 '13 at 10:14

1 Answer 1

up vote 1 down vote accepted

Could be interpreted as a nonlinear least squares problem.

pts = {{70.26, 45.78}, {71.04, 46.32}, {37.23, 24.67}, {37.91, 
28.78}};

ndsoln[{k1_, k2_, k3_}, x0_] := x[720] /. 
  NDSolve[{x'[t] == -k1 x[t]^2 + k2 x[t] y[t], 
  y'[t] == k1 x[t]^2 - k2 x[t] y[t] - k3 y[t], x[0] == x0, 
  y[0] == 0}, {x, y}, {t, 0, 720}][[1]]

ssfun[{k1_?NumericQ, k2_?NumericQ, k3_?NumericQ}] := 
  Sum[(ndsoln[{k1, k2, k3}, pts[[i, 1]]] - pts[[i, 2]])^2, 
  {i, Length[pts]}]

scaledArgs = {10^-5 k1, k2, 10^2 k3};

fm = FindMinimum[ssfun[scaledArgs], {k1, 3, 4}, {k2, 20, 21}, {k3, 5, 6}]

Show[Plot[ndsoln[scaledArgs /. fm[[2]], x0], {x0, 35, 75}], 
 ListPlot[pts, PlotStyle -> Directive[ColorData[1][2], PointSize[Medium]]], 
 PlotRange -> {20, 50}, Frame -> False, Axes -> True, 
 AxesLabel -> {"x(0)", "x(720)"}, AxesOrigin -> {35, 20}]

The result

share|improve this answer
    
Thank you very much! It was so helpful. –  Mai May 14 '13 at 0:43

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.