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