# How to Fit experimental data with 4D ParametricNDSolve?

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:

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

I have the differential equations with F(t)=mr″(t)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}\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}C_{drag} as a parameter; mm, \rho\rho, AA and RR 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*)

=================

Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
– Michael E2
May 6 ’15 at 23:45

– bbgodfrey
May 7 ’15 at 2:01

You are using Area in your code which is a Mathematica function
– Julian
May 7 ’15 at 8:25

Thanks for the comments, I changed Area to A (it wasn’t a problem) and added the values for the constants, hope you can help me now
– user3142983
May 7 ’15 at 10:09

Suggest you might find this a useful reference.
– Oleksandr R.
May 7 ’15 at 13:11

=================

2

=================

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!

Your code isn’t working for me… The first part is ok but the fit code gives me an error: The function value is not a number at {cdrag,Vo} = {1.91862,51.6635}.. I think the problem is in the NMinimize but I can’t understand how you manage to “have y as a function of x, and not t”
– user3142983
May 9 ’15 at 1:00

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 🙂

Why did you use Cdrag = 0.144135 ? I want to find the best value of Cdrag to fit the experimental data!
– user3142983
May 7 ’15 at 15:13

Only this (Cdrag = 0.13952;corrected ) value corresponds to the maximum distance Xmax=126.7 meters.In these equations, it is the only variable .Sorry my English is Google Translate.Maybe someone smarter solve your task 😉
– Mariusz Iwaniuk
May 7 ’15 at 15:28

I have a question?The ball is a tennis or golf or ? Cdrag is a constant depending on the velocity , roughness , viscosity ,Reynolds number.When the ball begins to rotate is created Mangus force,that this force has a major influence on the trajectory.
– Mariusz Iwaniuk
May 8 ’15 at 1:46

It’s a golf ball, Cdrag depends on the velocity and the velocity changes. Next I will implement the Magnus force, this is only the first part of the problem! Thanks for your contribution
– user3142983
May 9 ’15 at 0:50