1
$\begingroup$

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:

enter image description here

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*)
$\endgroup$
6
  • $\begingroup$ 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! $\endgroup$
    – Michael E2
    May 6, 2015 at 23:45
  • $\begingroup$ Please provide missing values for all your constants. $\endgroup$
    – bbgodfrey
    May 7, 2015 at 2:01
  • $\begingroup$ You are using Area in your code which is a Mathematica function $\endgroup$
    – Julian
    May 7, 2015 at 8:25
  • $\begingroup$ 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 $\endgroup$ May 7, 2015 at 10:09
  • $\begingroup$ Suggest you might find this a useful reference. $\endgroup$ May 7, 2015 at 13:11

2 Answers 2

3
$\begingroup$

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}]

Plot of model

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}]]

enter image description here

I hope this helps!

$\endgroup$
1
  • $\begingroup$ 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" $\endgroup$ May 9, 2015 at 1:00
1
$\begingroup$

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.

enter image description here

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.

enter image description here

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

Bye :)

$\endgroup$
5
  • $\begingroup$ Why did you use Cdrag = 0.144135 ? I want to find the best value of Cdrag to fit the experimental data! $\endgroup$ May 7, 2015 at 15:13
  • $\begingroup$ 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 ;-) $\endgroup$ May 7, 2015 at 15:28
  • $\begingroup$ 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. $\endgroup$ May 8, 2015 at 1:46
  • $\begingroup$ 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 $\endgroup$ May 9, 2015 at 0:50
  • $\begingroup$ Maybe this help GolfBall link, link $\endgroup$ May 9, 2015 at 11:43

Your Answer

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

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