2
$\begingroup$

I have a system of two ODEs depending on two parameter k,h and I want to find the best k,h to fit some data. To understand if my code works I started generating some fake data from the equations. Here is what I use to fit the data (I basically followed this post ODE fitting to dataset)

func=ParametricNDSolveValue[{x'[t]==k-x[t]/(1-x[t]^2-y[t]^2),y'[t]==h-y[t]/(1-x[t]^2-y[t]^2), x[0]==0,y[0]==0},{x,y},{t,0,3},{k,h}]
model[k_, h_][i_, t_] :=   Through[func[k, h][t], List][[i]] /; And @@ NumericQ /@ {k, h, i, t};
fitted = NonlinearModelFit[ dataset ,model[k, h][i, t],{{k,2/3}, {h,0.5}}, {i, t}] ;
Show[ListPlot[{Transpose[{time,datax}],Transpose[{time,datay}]}], Plot[ {fitted[1, t], fitted[2, t]}, {t, 0,3},PlotRange->All]]
fitted["ParameterTable"]

I generated my dataset solving the ODEs above with a specific value for k and h, and taking points only at some point in time. So I have the fake dataset that is in the form {index,time,data} (index 1 =data on x, index 2 =data on y) with time going from 0 to 3 with 0.1 steps.

{{1,0.,0.},{1,0.1,0.0316679},{1,0.2,0.0600199},{1,0.3,0.0851046},{1,0.4,0.107009},{1,0.5,0.12586},{1,0.6,0.14183},{1,0.7,0.155127},{1,0.8,0.165994},{1,0.9,0.174694},{1,1.,0.181503},{1,1.1,0.186694},{1,1.2,0.19053},{1,1.3,0.193256},{1,1.4,0.19509},{1,1.5,0.196221},{1,1.6,0.196812},{1,1.7,0.196998},{1,1.8,0.196888},{1,1.9,0.19657},{1,2.,0.196113},{1,2.1,0.19557},{1,2.2,0.19498},{1,2.3,0.194374},{1,2.4,0.193772},{1,2.5,0.193189},{1,2.6,0.192636},{1,2.7,0.192118},{1,2.8,0.191639},{1,2.9,0.191199},{1,3.,0.190799},{2,0.,0.},{2,0.1,0.0475773},{2,0.2,0.0905796},{2,0.3,0.129347},{2,0.4,0.164163},{2,0.5,0.195286},{2,0.6,0.222959},{2,0.7,0.247429},{2,0.8,0.268943},{2,0.9,0.287756},{2,1.,0.304121},{2,1.1,0.318291},{2,1.2,0.330508},{2,1.3,0.341005},{2,1.4,0.349996},{2,1.5,0.35768},{2,1.6,0.364233},{2,1.7,0.369814},{2,1.8,0.374562},{2,1.9,0.378599},{2,2.,0.38203},{2,2.1,0.384946},{2,2.2,0.387423},{2,2.3,0.389528},{2,2.4,0.391318},{2,2.5,0.39284},{2,2.6,0.394135},{2,2.7,0.395238},{2,2.8,0.396177},{2,2.9,0.396977},{2,3.,0.39766}} 

The code works but the fit is not completely perfected as I expected since I generated the data directly from equations (here is a picture).

Plot of data and fitted functioni.

I used k=1/3 and h=1/2 to generate data but it gives me as output k=0.26 and h =0.5018

If you want to run the code, here is the entire code:

f=1/3;
g=0.5;
sol = NDSolve[{x'[t]==f-x[t]/(1-x[t]^2-y[t]),y'[t]==g-y[t]/(1-x[t]^2-y[t]^2),x[0]==0,y[0]==0},{x,y},{t,0,3}]
datax=Flatten[Table[x[i]/.sol,{i,0,3,0.1}]];
datay=Flatten[Table[y[i]/.sol,{i,0,3,0.1}]];
time = Range[0,3,0.1];
datasetx=Transpose[{ConstantArray[1,31],time,datax}];
datasety=Transpose[{ConstantArray[2,31],time,datay}];
dataset=Join[datasetx,datasety]
func=ParametricNDSolveValue[{x'[t]==k-x[t]/(1-x[t]^2-y[t]^2),y'[t]==h-y[t]/(1-x[t]^2-y[t]^2), x[0]==0.0,y[0]==0.0},{x,y},{t,0,3},{k,h}]
model[k_, h_][i_, t_] :=   Through[func[k, h][t], List][[i]] /; And @@ NumericQ /@ {k, h, i, t};
fitted = NonlinearModelFit[ dataset ,model[k, h][i, t],{{k,0.3}, {h,0.4}}, {i, t}] ;
Show[ListPlot[{Transpose[{time,datax}],Transpose[{time,datay}]}], Plot[ {fitted[1, t], fitted[2, t]}, {t, 0,3},PlotRange->All]]
fitted["ParameterTable"]
$\endgroup$
4
  • $\begingroup$ Please provide your dataset too! $\endgroup$ Commented Dec 16, 2021 at 14:45
  • $\begingroup$ @UlrichNeumann I just edited the question adding more details on how i created the fake dataset. $\endgroup$ Commented Dec 16, 2021 at 14:52
  • $\begingroup$ I would expect dataset in the form {{t[i],{x[i],y[i]}},...}, but your code doesn't run. $\endgroup$ Commented Dec 16, 2021 at 15:00
  • $\begingroup$ As I said i followed what is done here mathematica.stackexchange.com/questions/84063/… which is what I need to do. I defined the model and the two index are used to refer to equation 1 or 2. I edited the question with the full code at the end that should run now. $\endgroup$ Commented Dec 16, 2021 at 15:11

1 Answer 1

2
$\begingroup$

With correct dataset your code & Mathematica evaluates as expected!

dataset = Join[Table[{1, t, func[1/3, 1/2][[1]][t]}, {t, 0, 3, .1}],
Table[{2, t, func[1/3, 1/2][[2]][t]}, {t, 0, 3, .1}]];

model[k_, h_][i_, t_] :=func[k, h][[i]][t] /; And @@ NumericQ /@ {k, h, i, t};
fitted = NonlinearModelFit[dataset,model[k, h][i, t], {{k, 2/3}, {h, 0.5}}, {i, t}];
fitted["BestFitParameters"]
(*{k -> 0.333333, h -> 0.5}*)
$\endgroup$
1
  • $\begingroup$ Yes, sorry, my bad! I just found a small error in the fake database. $\endgroup$ Commented Dec 17, 2021 at 11:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.