How to Fit experimental data with 4D ParametricNDSolve?

Question

Let me first explain the background problem. I'm trying to solve a problem with a projectile motion in 4D such that the position vector is:

$\boldsymbol{r}(t)=\boldsymbol{r}[x(t),y(t),z(t)]$

I have the differential equations with $F(t)=mr''(t)$ such that:

$\boldsymbol{F}_{drag}(t)=\frac{1}{2} \rho A C_{drag} R (\boldsymbol{r'}(t))^2; \\ \boldsymbol{F}_{grav}=m\boldsymbol{g};\\ \boldsymbol{F}_{total}(t)=\boldsymbol{F}_{drag}(t)+\boldsymbol{F}_{grav}$

I want to fit $C_{drag}$ as a parameter; $m$, $\rho$, $A$ and $R$ are known.

In Mathematica I'm solving them as:

r[t_] := {xx[t], yy[t], zz[t]};
(* Initial velocity V0 is known *)
V0 = {v0x,v0y,v0z};

(* Gravitical Force *)
GravForce = {0, -9.8*M, 0};
(* Drag Force *)
DragForce[t_] := 0.5*Cdrag*rho*A*Norm[r'[t]]*r'[t];
(* Total force *)
TotalForce[t_] := {DragForce[t][[1]] + GravForce[[1]], DragForce[t][[2]] + GravForce[[2]], DragForce[t][[3]] + GravForce[[3]]};

(* Parametric solution *)
ParametricSolution = ParametricNDSolve[{
r''[t][[1]] == TotalForce[t][[1]]/M, r'[0][[1]] == V0[[1]], r[0][[1]] == 0,
r''[t][[2]] == TotalForce[t][[2]]/M, r'[0][[2]] == V0[[2]], r[0][[2]] == 0,
r''[t][[3]] == TotalForce[t][[3]]/M, r'[0][[3]] == V0[[3]], r[0][[3]] == 0},
{xx, yy, zz}, {t, 20}, {Cdrag}];

Until here everything is working fine! But now I have a set of {x,y} experimental data that I want to find the best fit with Cdrag as a parameter.

data = {{0, 0}, {20, 8}, {40, 17}, {60, 24}, {80, 27}, {100, 
23}, {115, 14}, {126.7, 0}};

It looks like this:

I needed something like this:

NonlinearModelFit[data,
{xx[CMag][t], yy[CMag][t]} /. ParametricSolution,
{CMag}, {xx[t], yy[t]}];

But I get the error: "NonlinearModelFit::fitc: Number of coordinates (1) is not equal to the number of variables (2)."

Can you help me? Thanks ;)

EDIT:

The values for the constants are:

M = .04593; (* Mass of the ball in Kg *)
R = .04267/2; (* Ball Radius *)
A = Pi*R^2; (* Cross section area of the golf ball *)
rho = 1.2041; (* kg/m^3 Density of air at 20C, 1 atm*)
theta = 19.9;(* Starting angle from the horizontal *)
v0=48.54; (*speed in m/s*)
V0 = {v0*Cos[theta Degree], v0*Sin[theta Degree], 0};(* Ball velocity components*)

Answer 1

There is a problem with your initial velocity:

Hmax = Vo^2*Sin[phi Degree]^2/(2*9.8)
(*13.9274*)

With that initial Velocity, you can only reach 13.9 meters (far from the 25+ from your data), and that's without drag. So you can't fit anything really.

Maybe the velocity is not well measured. What you could do is fit Cdrag and Vo.

Using part of Mariusz's code

g = 9.80665;
M = 0.04593;
R = 0.04267/2;
phi = 19.9
rho = 1.2041;
A = Pi R^2;
k = 1/2 rho A cdrag;

sol = ParametricNDSolve[{x''[t] == -(k/M)*(x'[t])^2,
     y''[t] == -g - (k/M)*(y'[t])^2, x[0] == 0, 
    y[0] == 0, x'[0] == Vo*Cos[phi Degree], 
    y'[0] == Vo*Sin[phi Degree]}, {x, y}, {t, 0, 5},
    {cdrag, Vo}, Method -> "ExplicitRungeKutta"];

Then you can play around to see where are your parameters:

data = {{0, 0}, {20, 8}, {40, 17}, {60, 24}, {80, 27}, {100, 23},
        {115, 14}, {126.7, 0}};

sol2[t_, {cdrag_, Vo_}] := Evaluate[{x[cdrag, Vo][t], y[cdrag, Vo][t]} /. sol]

Manipulate[
   Show[
     ListPlot[data],
     ParametricPlot[sol2[t, {cdrag, Vo}], {t, 0, 5}]
   ], {cdrag, 0, 2}, {Vo, 30, 100}]

You can see your parameters are around 0.7 and 80 for cdrag and Vo

Now, to produce a fit, I will do the following:

(*To have y as a function of x, and not t*)
yval[xval_, {cdrag_, Vo_}] := Module[
  {tval},
  tval = t /. FindRoot[Evaluate[x[cdrag, Vo][t] /. sol] == xval, {t, 3}] //Quiet;
  Evaluate[y[cdrag, Vo][tval] /. sol]
  ]

(*This is to get the y values of the model at the x values of data*)
modelPoints[cdrag_?NumericQ, Vo_?NumericQ] := 
  Map[{#, yval[#, {cdrag, Vo}]} &, data[[ ;; , 1]]]

(*Using NMinimize to fit by Least Squares*)
fit = NMinimize[
   {(modelPoints[cdrag, Vo] - data)^2 // Total // Last,
    0 < cdrag < 2, Vo > 50}
   , {cdrag, Vo}] // Quiet

(*{cdrag -> 0.70642, Vo -> 79.7468}*)


{fittedCdrag, fittedVo} = fit[[All, 2]];
Show[
 ListPlot[data],
 ParametricPlot[sol2[t, {fittedCdrag, fittedVo}], {t, 0, 5}]]

I hope this helps!

Answer 2

OK, I'm able to give an incomplete answer now.I can't fit but..

Values for all constants are correct? You must have made ​​a mistake somewhere.

Cdrag = 0.13952;
g = 9.80665;
M = 0.04593;
R = 0.04267/2;
Vo = 48.54;
\[Phi] = 19.9;
\[Rho] = 1.2041;
A = \[Pi] R^2;
k = 1/2 \[Rho] A Cdrag;

{xsol, ysol} = NDSolveValue[{x''[t] == -(k/M)*(x'[t])^2, 
y''[t] == -g - (k/M)*(y'[t])^2, x[0] == 0, y[0] == 0, 
x'[0] == Vo*Cos[\[Phi] Degree], 
y'[0] == Vo*Sin[\[Phi] Degree]}, {x, y}, {t, 0, 5}, 
Method -> "ExplicitRungeKutta"];

{s = FindRoot[ysol[t] == 0, {t, 5}], {"X max =", xsol[t] /. s}}

Output:{{t -> 3.29119}, {"X max =", 126.7}} with yours data y=0 is Xmax=126.7 meters.

{Hmax = FindRoot[ysol'[t] == 0, {t, 5}], ysol[t] /. Hmax}

Output:{{t -> 1.64559}, 13.4347} Hmax is=13,4347 meters

Plot.

Show[ListPlot[data, PlotStyle -> Black, AxesLabel -> {X[meters],Y[meters]}], ParametricPlot[{xsol[t], ysol[t]}, {t, 0, 4}, PlotStyle -> Red]]

Bye :)