The goal of the present problem is to be able to estimate the parameters $k$, $b$ and $m$ based on lists of known boundary conditions.
This is what I have so far
(*Constants*)
g = 9.81;
theta = {20,30,40}*(Pi/180);
phi = {45,90,135}*(Pi/180);
v = 35;
Ux = {5,2,8};
Uy = {0,2,3};
Uz = ConstantArray[0,3];
(*At `z==0` (coordinates in (x,y,z))*)
P = {{10, 25, 0},{15,20,0},{16,18,0}};
(*Solve ODEs*)
S = Integrate[Sqrt[x'[t]^2 + y'[t]^2 + z'[t]^2], t];
sol[v_, theta_, phi_, {Ux_, Uy_, Uz_}] =
ParametricNDSolve[{z''[t] == -m*g -
Tanh[z'[t]]*kz (z'[t] - Uz)^2 (1 + Exp[bz*S]), z[0] == 4,
z'[0] == v*Cos[theta],
x''[t] == -Tanh[x'[t]]*kx (x'[t] - Ux)^2 (1 + Exp[bx*S]),
x[0] == 0, x'[0] == v*Sin[theta]*Cos[phi],
y''[t] == -Tanh[y'[t]]*ky (y'[t] - Uy)^2 (1 + Exp[by*S]),
y[0] == 0, y'[0] == v*Sin[phi]*Sin[theta],
WhenEvent[z[t] == 0, {tMax = t, "StopIntegration"}]}, {x, y,
z}, {t, 0, EndTime[theta]}, {m, kx, ky, kz, bx, by, bz},
Method -> "StiffnessSwitching"];
The next step is to estimate the parameters such that boundary conditions are met as accurately as possible. What is the best way to make the solutions of ParametricNDSolve fit lists of BCs?
Using NMinimizeworks for a single configuration
(*Constants*)
g = 9.81;
theta = 45*(Pi/180);
phi = 45*(Pi/180);
v = 35;
Ux = 5;
Uy = 5;
Uz = 0;
(*At `z==0` (coordinates in (x,y,z))*)
P = {70, 70, 0};
(*Solve ODEs*)
S = Integrate[Sqrt[x'[t]^2 + y'[t]^2 + z'[t]^2], t];
sol[v_, theta_, phi_, {Ux_, Uy_, Uz_}] =
ParametricNDSolve[{z''[t] == -m*g -
Tanh[z'[t]]*kz (z'[t] - Uz)^2 (1 + Exp[bz*S]), z[0] == 0,
z'[0] == v*Cos[theta],
x''[t] == -Tanh[x'[t]]*kx (x'[t] - Ux)^2 (1 + Exp[bx*S]),
x[0] == 0, x'[0] == v*Sin[theta]*Cos[phi],
y''[t] == -Tanh[y'[t]]*ky (y'[t] - Uy)^2 (1 + Exp[by*S]),
y[0] == 0, y'[0] == v*Sin[phi]*Sin[theta],
WhenEvent[z[t] == 0, {tMax = t, "StopIntegration"}]}, {x, y,
z}, {t, 0, EndTime[theta]}, {m, kx, ky, kz, bx, by, bz},
Method -> "StiffnessSwitching"];
Quiet@Minimize[
EuclideanDistance[P,
Through[#[#[[1, 1, -1, -1]]]] &@
Through[({x, y, z} /. sol[v, theta, phi, {Ux, Uy, Uz}])[m, kx, ky,
kz, bx, by, bz]]], {m, kx, ky, kz, bx, by, bz}]
(*{0.00013925, {m -> 0.121799, kx -> 0.139893, ky -> 0.0000576328,
kz -> 0.267708, bx -> -0.465209, by -> 0.406584, bz -> 0.0751821}}*)
Background
The above ODEs describes (empirically) the motion of a liquid jet in the presence of air drag(1). The drag force is
$$kv^2(1+\exp (bS)$$
and $S$ is the total travelled distance by the jet
$$S=\int\sqrt{x'[t]^2+y'[t]^2+z'[t]^2}dt$$
The reason for the exponential term is that initially upon exit from the nozzle, a water jet is seemingly unaffected by air drag, but quickly disintegrates over long distances after it has been destabilized.
$P$ is measurements of landing points of a liquid jet at configurations corresponding to angles $\varphi$ and $\theta$ (spherical coordinates) under the influence of air with speed $U$.
