3
$\begingroup$

Cross-posted in http://tieba.baidu.com/p/4440556337 http://tieba.baidu.com/p/4440547130

I was wondering if there is a function that fits the following data perfectly or properly. There is a requirement that the function must be a decreasing function when the independent variable t is in range 0 to 5. The data is as follows:

{{0, 282.843}, {1, 18.071}, {2, 6.2538}, {3, -3.58047}, {4, -12.6429}, {5, -18.3376}}

Ok,the source of my data is the numerial solution of the differential equation of vy1,and the codes are as follows:

Remove["Global`*"];
θ = π/4; g = 9.8; v0 = 400; k = 0.1; m = 100;
v0x = v0*Cos[θ]; v0y = v0*Sin[θ];
vx1 = DSolve[{-k*(vx[t])^3 == m*vx'[t], vx[0] == v0x}, vx[t], t]
vy1 = NDSolve[{-k*(vy[t])^3 - m*g == m*vy'[t], vy[0] == v0y}, 
  vy, {t, 0, 5}]
Plot[{vy[t] /. vy1}, {t, 0, 5}]
data = Table[Flatten[{t, vy[t]} /. vy1], {t, 0, 5}]

and you can get the shape of the curve of vy1 from the plot,so how can we get the best fitting of vy1,thanks!

$\endgroup$
  • 1
    $\begingroup$ This really is more a math/statistics question than a Mathematica one. Nevertheless: where did the data come from? Is there no theory behind that can supply a possible model? $\endgroup$ – J. M. will be back soon Mar 27 '16 at 12:34
  • $\begingroup$ @J.M. I have just added my supplement and may it is useful to help solve the question,thanks! $\endgroup$ – dcydhb Mar 27 '16 at 13:48
7
$\begingroup$

Update

Updated to handle the cases k=0.1, k=1 and k=10.

Using your parameters

θ = π/4;
g = 9.8;
v0 = 400;
m = 100;
v0y = v0*Sin[θ];

k=0.1

k=0.1;
vy1kp1 = NDSolve[{-k*(vy[t])^3 - m*g == m*vy'[t], vy[0] == v0y}, vy, {t, 0, 5}]
vykp1 = vy1kp1[[1, 1, 2]];

k=1

k = 1;
vy1k1 = NDSolve[{-k*(vy[t])^3 - m*g == m*vy'[t], vy[0] == v0y}, vy, {t, 0, 5}]
vyk1 = vy1k1[[1, 1, 2]];

k=10

k = 10;
vy1k10 = NDSolve[{-k*(vy[t])^3 - m*g == m*vy'[t], vy[0] == v0y}, vy, {t, 0, 5}]
vyk10 = vy1k10[[1, 1, 2]];

Plot the data

Plot[{vykp1[t], vyk1[t], vyk10[t]}, {t, 0, 5},
 PlotRange -> {{0, 5}, {-20, 69}}, PlotStyle -> {Blue, Black, Red}]

Mathematica graphics

After analyzing the data it appears that the early part could be modeled as two exponential declines (one fast and the second slower) followed by a linear portion. The late data asymtotes to a constant.

We need a function to blend between the early and late part

blendFun[t_, tblend_, texp_, dt_: 10^-8] :=
 (1 + ((tblend/(t + dt))^texp - ((t + dt)/tblend)^
        texp)/((tblend/(t + dt))^texp + ((t + dt)/tblend)^texp))/2


Plot[{blendFun[t, 3, 5], (1 - blendFun[t, 3, 5])}, {t, 0, 5}, 
 PlotStyle -> {Black, Blue}, PlotRange -> {{0, 5}, {0, 1}}]

Mathematica graphics

k=1 case

This is the simplest case in the sense that there is sufficient and approximately equal amounts of early and late time data that makes the fitting operation straight forward.

Constraints are imposed to set the sign of the parameters. In addition a constraint is set to force the sum of the two exponentials and the linear offset to equate to the first data value. Finally the constant term is forced to be equal to the final data value.

Clear[m]
datak1 = Table[Flatten[{t, vyk1[t]}], {t, 0, 5, 0.05}];

solk1 = With[
  {
   blend = ((1 + ((tblend/(t + 10^-6))^texp - ((t + 10^-6)/tblend)^
            texp)/((tblend/(t + 10^-6))^texp + ((t + 10^-6)/tblend)^
            texp))/2)
   },

  FindFit[
   datak1, {a Exp[-b t] + c Exp[-d t] + (o + m t) blend + 
     constant (1 - blend), {a > 0, b > 0, c > 0, d > 0, o > 0, m < 0, 
     tblend > 0, texp > 1, constant == datak1[[-1, 2]], 
     a + c + o == datak1[[1, 2]]}}, {{a, 250}, {b, 75.0}, {c, 20}, {d,
      5}, {m, -5.0}, {o, 9.}, {tblend, 2.0}, {texp, 
     2}, {constant, datak1[[-1, 2]]}}, t, MaxIterations -> 1000]
  ]

{a -> 252.335, b -> 74.9384, c -> 21.8433, d -> 5.12808, 
 m -> -5.33564, o -> 8.66416, tblend -> 1.78099, texp -> 2.12587, 
 constant -> -9.93229}

Show[ListLinePlot[datak1, PlotStyle -> Black],
 Plot[
  Evaluate[
   With[
    {
     blend = ((1 + ((tblend/(t + 10^-6))^texp - ((t + 10^-6)/tblend)^
              texp)/((tblend/(t + 10^-6))^texp + ((t + 10^-6)/tblend)^
              texp))/2)
     },
    a Exp[-b t] + c Exp[-d t] + (o + m t) blend + 
      constant (1 - blend) /. solk1]
   ], {t, 0, 5}, PlotStyle -> Red
  ],
 PlotRange -> All
 ]

Mathematica graphics

k = 0.1 case

This is more difficult as there is not sufficient late time data to locate the constant value. We will impose a constraint that forces the constant to be less than the last data value. That is the best that we can do in this case.

datakp1 = Table[Flatten[{t, vykp1[t]}], {t, 0, 5, 0.1}];

solkp1 = With[
  {
   blend = ((1 + ((tblend/(t + 10^-6))^texp - ((t + 10^-6)/tblend)^
            texp)/((tblend/(t + 10^-6))^texp + ((t + 10^-6)/tblend)^
            texp))/2)
   }, 

  FindFit[
   datakp1, {a Exp[-b t] + c Exp[-d t] + (o + m t) blend + 
     constant (1 - blend), {a > 0, b > 0, c > 0, d > 0, o > 0, m < 0, 
     tblend > 0, texp > 1, constant <= datakp1[[-1, 2]], 
     a + c + o == datakp1[[1, 2]]}}, {{a, 210}, {b, 30.0}, {c, 
     50}, {d, 4}, {m, -9.0}, {o, 25.}, {tblend, 5.0}, {texp, 
     2}, {constant, -20}}, t, MaxIterations -> 1000]
  ]

{a -> 210.918, b -> 30.022, c -> 46.7299, d -> 3.34615, m -> -8.80741,
  o -> 25.1946, tblend -> 5.2677, texp -> 1.55721, 
 constant -> -18.3376}

Show[ListLinePlot[datakp1, PlotStyle -> Black], 
 Plot[a Exp[-b t] + c Exp[-d t] + m t + o + e Exp[s t] /. solkp1, {t, 
   0, 5}, PlotStyle -> Red], PlotRange -> All]

Mathematica graphics

k = 10 case

This is also difficult. There is much more late time than early time data so the unbalance makes it difficult to fit in one fell swoop.

Instead we will fit the early time data first and then impose the exponential terms when fitting the late time data.

datak10 = Table[Flatten[{t, vyk10[t]}], {t, 0, 5, 0.001}];
solk10 = FindFit[
  datak10[[1 ;; 
    701]], +{a Exp[-b t] + c Exp[-d t] + o + m t, {a > 0, 
     2000 > b > 0, c > 0, d > 0, o > 0, -20 < m < 0, 
     a + c + o == datak10[[1, 2]]}}, {{a, 245}, {b, 1900.0}, {c, 
    30}, {d, 60}, {m, -12}, {o, 7.}}, t, MaxIterations -> 1000]

{a -> 244.029, b -> 1957.1, c -> 31.6787, d -> 57.0243, m -> -12.8196,
  o -> 7.1351}

We will impose the results for the exponential components as well as the constant term when fitting the linear portion and the blend. This can be done as follows:

solk10 = With[
  {
   a = 244.02887,
   b = 1957.1,
   c = 31.6787,
   d = 57.0243,
   constant = datak10[[-1, 2]],
   blend = ((1 + ((tblend/(t + 10^-6))^texp - ((t + 10^-6)/tblend)^
            texp)/((tblend/(t + 10^-6))^texp + ((t + 10^-6)/tblend)^
            texp))/2)
   },

  FindFit[
   datak10, {a Exp[-b t] + c Exp[-d t] + (o + m t) blend + 
     constant (1 - blend), {-15 < m < -5, o > 0, tblend > 0.5, 
     2 > texp > 1}}, {{m, -6}, {o, 7}, {tblend, 0.6}, {texp, 1.2}}, t,
    MaxIterations -> 1000]
  ]

{m -> -5.85949, o -> 6.71999, tblend -> 0.588783, texp -> 1.20832}

And now plot the fit

Show[ListLinePlot[datak10, PlotStyle -> Black, PlotRange -> {{0, 5}, {-6, 39}}],
 Plot[
  Evaluate[
   With[
    {
     a = 244.02887,
     b = 1957.1,
     c = 31.6787,
     d = 57.0243,
     (*m=-12.8196,
     o=7.1351,*)
     constant = datak10[[-1, 2]],
     blend = ((1 + ((tblend/(t + 10^-6))^texp - ((t + 10^-6)/tblend)^
              texp)/((tblend/(t + 10^-6))^texp + ((t + 10^-6)/tblend)^
              texp))/2)
     },
    a Exp[-b t] + c Exp[-d t] + (o + m t) blend + 
      constant (1 - blend) /. solk10]
   ], {t, 0, 5}, PlotStyle -> Red, PlotRange -> {{0, 5}, {-6, 39}}
  ]
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Always excited by the wonderful answers and also this time,surely a great answer but there is also some points that I was curious about,for example how to set the initial values from a to s, ` {{a, 212}, {b, 30.0}, {c, 50}, {d, 3}, {m, -10.0}, {o, 25.}, {e, 0.01}, {s, 1.1}} ,for example in this case vy2 = NDSolve[{-k*(vy[t])^3 + mg == mvy'[t], vy[0] == 0}, vy, {t, 0, 10}] Plot[{vy[t] /. vy2}, {t, 0, 10}] data2 = Table[Flatten[{t, vy[t]} /. vy2], {t, 0, 10, 0.1}];`I was wondering how to set the initial values ,thanks a lot! $\endgroup$ – dcydhb Mar 28 '16 at 8:19
  • $\begingroup$ And also in my differential equation the value of ` k = 0.1,but when k=1,k=10`,then what is the initial value from a to s,thanks! $\endgroup$ – dcydhb Mar 28 '16 at 9:07
  • $\begingroup$ I found the initial values by wrapping the Show in a Manipulate and manually adjusting the sliders until I found an approximate fit. Manipulate[ Show[ ListLinePlot[data, PlotStyle -> Black], Plot[a Exp[-b t] + c Exp[-d t] + m t + o , {t, 0, 5}, PlotStyle -> Red], PlotRange -> All ], {{a , 200}, 10, 400, Appearance -> "Open"}, {{b , 30}, 15, 40, Appearance -> "Open"}, {{c, 50}, 10, 100, Appearance -> "Open"}, {{d, 3}, 1, 5, Appearance -> "Open"}, {{m, -10}, -30, 0, Appearance -> "Open"}, {{o , 25}, 0, 60, Appearance -> "Open"} ] $\endgroup$ – Jack LaVigne Mar 28 '16 at 13:09
  • $\begingroup$ @ Jack LaVigne I have tried whenk=10 and {{a, 231.5}, {b, 83.5}, {c, 62.1}, {d, 0.674}, {n, 17}, {o, -53}, {e, 0.01}, {s, 1.1}},but there is still some deviation existing and not a perfect fitting. $\endgroup$ – dcydhb Mar 28 '16 at 14:38
  • $\begingroup$ To cover the cases where k is set to 0.1, 1 and 10 a different form is needed. See the update. $\endgroup$ – Jack LaVigne Apr 2 '16 at 15:02
4
$\begingroup$

If you have Mathematica 10.2 or higher, you can try the new FindFormula function, which uses symbolic regression to search for formulae fitting your data using target functions you think would be likely to appear in a solution.

FindFormula[{{0, 282.843}, {1, 18.071},...},t,TargetFunctions->{Plus,Times,Log,Sin,...}]

It might also help to generate more than five data points - If you generate the data with data = Table[Flatten[{t, vy[t]} /. vy1], {t, 0, 5, n}], you can set n to be a much smaller increment than 1 and get far more points, which you'll need regardless of what method or program you use to try to fit them.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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