2

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"]
user_1922487
  • 23
  • 3
  • Please provide your dataset too! – Ulrich Neumann Dec 16 '21 at 14:45
  • @UlrichNeumann I just edited the question adding more details on how i created the fake dataset. – user_1922487 Dec 16 '21 at 14:52
  • I would expect dataset in the form {{t[i],{x[i],y[i]}},...}, but your code doesn't run. – Ulrich Neumann Dec 16 '21 at 15:00
  • As I said i followed what is done here https://mathematica.stackexchange.com/questions/84063/ode-fitting-to-dataset?newreg=2dcff3e4ffb942adb2ebad4f5aaa9006 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. – user_1922487 Dec 16 '21 at 15:11

1 Answers1

2

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})
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55