Model calibration with phase space data

Question

Forgive me if this is answered elsewhere, I am new to Mathematica and this forum and have tried a number of search phrases to no avail. My question is very general but for example, suppose I throw a ball into the air and record sequential positions on the same image (a motion diagram) but my camera faulted and didn't open the shutter in equal intervals. I can only measure y(x) rather than x(t) and y(t).

In Mathematica I can assume a physical model (quadratic drag) and, using ParametricNDSolve, I find numerical solutions to x(t) and y(t) with initial speed, v0, and drag, d, as fitting parameters.

I can easily fit x(t) or y(t) to the respective data but I am restricted to fitting y(x) to my available data. Armed with numerical solutions to x(t) and y(t), how would I find a numerical fit to my y(x) data? I can plot y(x) just fine but how do I fit the data? My example code resulting in x(t) and y(t) with fabricated data is below.

Clear[d, v0]; g = 9.81;ϕ0 = 40 Degree;
d::usage = 
  "d has units of length and is the ratio of inertia over viscosity";
v0::usage = "v0 is the initial speed with upwards considered positive";
ϕ0::usage = "ϕ0 is the angle at which the ball is thrown";
eqnx = {x''[t] + (1/d)*Sqrt[(x'[t])^2 + (y'[t])^2]*x'[t] == 0};
eqny = {y''[t] + g + (1/d)*Sqrt[(x'[t])^2 + (y'[t])^2]*y'[t] == 0};
ic = {x'[0] == v0*Cos[ϕ0], x[0] == 0, y'[0] == v0*Sin[ϕ0], y[0] == 2};
solPN = ParametricNDSolve[{eqnx, eqny, ic}, {x, y}, {t, 0, 6}, {d, v0}]
data = Table[{t, (x[8, 23][t] /. solPN) + 
   RandomInteger[{-5, 5}]/30, (y[8, 23][t] /. solPN) + 
   RandomInteger[{-5, 5}]/30}, {t, 0.1, 3, 0.2}];

Thank you for your help.

Answer 1

There are several ways to do this. I am going to show a more-or-less straightforward way. Some adjustments probably have to be done when using real data.

General procedure

1. Using the data find an interpolation (or use a fit with appropriate basis) of y[x].

2. Define a function ParDist that measures the deviation of the fit. That function takes as arguments:

   • the found parametric functions,

   • a pair of parameters,

   • an interpolation function for y[x],

   • a vector of time points (to be used to compute the measure).

3. Select a set of time points for ParDist.

4. Minimize the measure function ParDist over an appropriate domain for the parameters.

5. Compute new data points with the found parameters over the selected time points.

6. Compare the new data with the initial data.

7. Repeat the steps 1-6 using different minimization methods, norms in ParDist, and sets of time points for ParDist.

It might happen that the interpolation in step 1 is not an adequate approach, and some fitting procedure has to be used. For simplicity I am using interpolation below.

Code

Interpolation function for y[x]:

xyFunc = Interpolation[data[[All, 2 ;; 3]]]

The measure function to be minimized, ParDist:

Clear[ParDist]
ParDist[{x_ParametricFunction, y_ParametricFunction}, {d_?NumberQ, v0_?NumberQ}, xyFunc_, ts : {_?NumberQ ..}] :=
  Block[{points, xMin, xMax},
   points = Map[{x[d, v0][#], y[d, v0][#]} &, ts];
   {xMin, xMax} = xyFunc["Domain"][[1]];
   points = Select[points, xMin <= #[[1]] <= xMax &];
   If[Length[points] == 0, 10^10.,
    Norm[Map[(xyFunc[#[[1]]] - #[[2]])/Abs[xyFunc[#[[1]]]] &, points]]
   ]
  ];

(The selection of the points to be within the interpolation domain of xyFunc is important.)

Select a set of time points (does not have to be a regular grid):

ts = Range[0.5, 3, 0.13];

Minimize ParDist over an appropriate domain for the parameters:

nsol = 
 NMinimize[{ParDist[{x, y} /. solPN, {d, v0}, xyFunc, ts], 0 <= d <= 20, 5 <= v0 <= 25}, {d, v0}]

(* {0.654588, {d -> 9.2491, v0 -> 19.9544}} *)

Compute new data points with the found parameters:

Block[{d, v0},
  {d, v0} = {d, v0} /. nsol[[2]];
  fitData = Table[{t, x[d, v0][t] /. solPN, y[d, v0][t] /. solPN}, {t, ts}]
];

Compare the plots of the fitted data and initial data:

ListPlot[{data[[All, 2 ;; 3]], fitData[[All, 2 ;; 3]]}, 
 PlotRange -> All, PlotTheme -> "Detailed", 
 PlotLegends -> {"initial", "fitted"}]

Answer 2

You can take the inverse of x[t] and feed it back into y[t] to get a y[x] directly. We can do this because for this particular problem formulation $x^{-1}(t)$ is single valued. If we were shooting into a headwind where the projectile might eventually curve back around, all bets would be off.

Use @jfaraday's code but replace the ParametricNDSolve[ ] statement.

{xtt, ytt} = ParametricNDSolve[{eqnx, eqny, ic}, {x, y}, {t, 0, 20}, {d, v0}] // Values

Note the time runs out to 20 seconds, to avoid having to extrapolate during the optimization.

In general we can express y[x] as

yx[xInput_]:=ytt[dd, v00][InverseFunction[xtt[dd, v00]][xInput]]

Now solve, but first splitting out some data to make things clear. Also, using FindMinimum[ ] with reasonable start values, stolen from @Anton.

xData = data[[All, 2]];
yData = data[[All, 3]];

res = FindMinimum[Norm[yData - (ytt[dd, v00][InverseFunction[xtt[dd, v00]][#]] & /@ xData)], 
                 {{dd, 9},{v00, 20}}]

(* {1.25097, {dd -> 8.24338, v00 -> 22.3428}} *)

Plot it

yfit = (ytt[dd, v00][InverseFunction[xtt[dd, v00]][#]] & /@ xData) /. res[[2]]

ListPlot[{data[[All, {2, 3}]], Transpose@{xData, yfit}}]