In[]:=
g=–9.81;v0=0;h0=100000;k=0.001;m0=2000;T=15000;
In[]:=
pfun=ParametricNDSolveValue[{z''[t]g+T/m[t],z'[0]v0+gτ,z[0]h0+v0τ+(1/2)gτ^2–s[t]^2,m'[t]–kT,m[0]m0,s[0]0},z,{t,0,120},{τ}]
Out[]=
ParametricFunction
 |  |
|
In[]:=
fullZ[τ_,t_,s_]:=If[t>τ,pfun[τ][t–τ,s],h0+v0t+(1/2)gt^2]fullV[τ_,t_,s_]:=If[t>τ,pfun[τ]'[t–τ,s],v0+gt]
In[]:=
Plot[Evaluate[Table[fullZ[τ,t,0],{τ,70,80,2}]],{t,0,200},PlotRangeAll]
Out[]=
root=FindRoot[{fullZ[τ,tl,s],fullV[τ,tl,s]},{{τ,75},{tl,190},{s,0}}]
Out[]=
{τ75.7512,tl187.999}
In[]:=
Show[Plot[Evaluate@Table[Tooltip[fullZ[τ,t],τ],{τ,70,80,2}],{t,0,200},PlotStyleLighter[Gray,0.5]],Plot[Evaluate[fullZ[τ,t]/.root],{t,0,Evaluate[tl/.root]},PlotStyleOrange]]
Out[]=
In[]:=
TableForm[Table[Evaluate[{t,fullZ[τ,t],fullV[τ,t]}/.root],{t,183,193.,0.5}],TableHeadings{None,{"t","h(t)","v(t)"}}]
Out[]//TableForm=
t | h(t) | v(t) |
183. | 428.076 | –163.721 |
183.5 | 349.812 | –149.271 |
184. | 278.868 | –134.439 |
184.5 | 215.439 | –119.211 |
185. | 159.727 | –103.568 |
185.5 | 111.943 | –87.4942 |
186. | 72.3074 | –70.9709 |
186.5 | 41.0501 | –53.9784 |
187. | 18.4107 | –36.4956 |
187.5 | 4.63991 | –18.5001 |
188. | 0.000013704 | 0.0321097 |
188.5 | 4.76576 | 19.1269 |
189. | 19.2252 | 38.8118 |
189.5 | 43.6809 | 59.1167 |
190. | 78.4506 | 80.0736 |
190.5 | 123.869 | 101.717 |
191. | 180.288 | 124.085 |
191.5 | 248.081 | 147.217 |
192. | 327.64 | 171.158 |
192.5 | 419.382 | 195.956 |
193. | 523.748 | 221.664 |
Leave a Reply