Message Boards Message Boards

GROUPS:

Estimating parameters based on lists of boundary conditions

Posted 8 years ago
5219 Views
|
1 Reply
|
0 Total Likes
|

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$.

POSTED BY: Julian Lovlie
Posted 8 years ago

I understand the goal but given the code listed, haven't you already achieved it? For any set of measurements and landing points you can already calculate/estimate k and m.

If you you have multiple sets measurements and landing points that result from the same "true" values of k and m, then you could apply your code as a function, obtain the multiple estimates of k and m and even calculate measures of precision. (With the measures of precision associated with how you obtain multiple sets of measurements and landing points.)

If the multiple sets of measurements and landing points are from different "true" values of k and m, then you could try to estimate k and m using a simpler function of the measurements (using NonlinearModelFit). But you also need to specify the error structure of the measurements and landing points - i.e., how the experimental data is generated. Is there a single projectile? Or multiple projectiles with different masses and/or drag constants? How is the start speed varied among runs of the experiment?

And, finally (almost done with the sermon), the variability of the estimates might very well differ depending on the projectile's properties (lighter projectiles might have more variability than heavier projectiles, for example). It's not clear to me that NonlinearFit allows much flexibility in setting the covariance structure.

POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract