3

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

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Mai
  • 33
  • 2
  • ParametricNDSolve[] is what seems to be used now, if you have version 9. – J. M.'s missing motivation May 11 '13 at 02: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 06: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 Answers1

2

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

mef
  • 1,629
  • 11
  • 15