3
$\begingroup$

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

The mathematica helps me solve one differential equation in numerial method and got the data of the answer,but the taylor series fit is not good in this situation, so I try to use FindFit to get the best concrete expression,but I was wondering why I can't get the result?Thanks!

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, 30}]
Plot[{vy[t] /. vy1}, {t, 0, 30}]
data = Table[Flatten[{t, vy[t]} /. vy1], {t, 0, 30}]
vy1' = Fit[ data, {1, t, t^2, t^3, t^4, t^5, t^6, t^7, t^8, t^9, t^10}, t]
Plot[{ vy1'}, {t, 0, 30}, Epilog -> Map[Point, data]]

model = a*Log[t] + b + c*Exp[t];
fit = FindFit[data, model, {a, b, c}, t]
Plot[Evaluate[model /. fit], {t, 0, 30}]

@bbgodfrey Certainly also a good supplement answer this time and I have tested and vertified that wheather k=0.1,k=1or k=10,it's accurate and almost no deviation ,but excuse me once again there are also some questions that puzzle me such as FIRST when we use taylor series fitting and we can get the plot of a large scale while in your method e can only get the curve in small scale ,so can you change it to a larger scale , NEXT, in your method when k=0.1,I got the concere formula expression of vy1,and it is in InverseFunction expression,as we all know that the integral of vis y ,but when I Integrate vy1 from 0 to a variable such as tau,the we can't get the concere expression or the formula of y with which I can get the plot of y-t,so how can we solve this problem THIRD,if we can get a correct or proper model to fit the expression or formula of vy1,for example using findfit or any other method,thus we can get the concrete expression or formula of y thanks a lot!

$\endgroup$
8
$\begingroup$

Solution for k = .1

The code in the question evaluates vy[t] numerically with NDSolve and then attempts to obtain an analytical approximation to it. The obvious alternative of using DSolve instead is not successful.

DSolve[{-k*(vy[t])^3 - m*g == m*vy'[t], vy[0] == v0y}, vy, {t, 0, 30}]
(* DSolve::bvnul: For some branches of the general solution, the given 
       boundary conditions lead to an empty solution. >> 
   {} *)

However, the error message hints at using DSolve without a boundary condition.

vy /. DSolve[{-k*(vy[t])^3 - m*g == m*vy'[t]}, vy, t][[1, 1]] /. C[1] -> c
(* Function[{t}, InverseFunction[(2 Sqrt[3] ArcTan[(-35 + 35^(1/3) #1)/(35 Sqrt[3])] + 
   2 Log[70 + 35^(1/3) #1] - Log[4900 - 70 35^(1/3) #1 + 35^(2/3) #1^2])/
   (840 35^(1/3)) &][-0.001 t + c]] *)

Next, obtain the value of c corresponding to vy[0] == v0y

FindRoot[%[0] == v0y, {c, .0016, 0, .00198}]
(* {c -> 0.0019741} *)

and insert it into the Function above for vy[t].

fvy = %% /. %

Not surprisingly, it agrees well with data, computed in the question.

Show[Plot[fvy[t], {t, 0, 30}, PlotRange -> {Automatic, {40, -22}}], 
    AxesLabel -> {t, vy}], ListPlot[data]]

enter image description here

This accurate solution, fvy[t], can be used like any other function of t.

Addendum: Solution for general k, g, m

At the request of dcydhb in a comment below, here is a more general approach.

Clear[k, m, g];
s = vy /. DSolve[{-k*(vy[t])^3 - m*g == m*vy'[t]}, vy, t][[1, 1]]
(* Function[{t}, InverseFunction[(
   2 Sqrt[3] ArcTan[(-1 + (2 k^(1/3) #1)/(g^(1/3) m^(1/3)))/Sqrt[3]] + 
   2 Log[g^(1/3) m^(1/3) + k^(1/3) #1] - Log[g^(2/3) m^(2/3) - g^(1/3) k^(1/3) m^(1/3) #1
   + k^(2/3) #1^2])/(6 g^(2/3) k^(1/3) m^(2/3)) &][-(t/m) + C[1]]] *)

For t -> 0, the inverse of this InverseFunction applied to v0y is C[1].

cs = s[[2, 0, 1]]
(* (2 Sqrt[3] ArcTan[(-1 + (2 k^(1/3) #1)/(g^(1/3) m^(1/3)))/Sqrt[3]] + 
   2 Log[g^(1/3) m^(1/3) + k^(1/3) #1] - Log[g^(2/3) m^(2/3) - 
   g^(1/3) k^(1/3) m^(1/3) #1 + k^(2/3) #1^2])/(6 g^(2/3) k^(1/3) m^(2/3)) & *)

For instance,

cs[v0y] /. {g -> 9.8, k -> 1/10, m -> 100}
(* 0.019741 *)
cs[v0y] /. {g -> 9.8, k -> 10, m -> 100}
(* 0.00426591 *)

Then,

Show[Plot[Evaluate[s[t] /. (C[1] -> cs[v0y]) /. {g -> 9.8, k -> 1/10, m -> 100}], 
    {t, 0, 30}, PlotRange -> {Automatic, {40, -22}}, AxesLabel -> {t, vy}], ListPlot[data]]

reproduces the Plot above, and

Plot[Evaluate[s[t] /. (C[1] -> cs[v0y]) /. {g -> 9.8, k -> 10, m -> 100}], 
    {t, 0, 4}, PlotRange -> {Automatic, {15, -5}}, AxesLabel -> {t, vy}]

yields the solution for k = 10.

enter image description here

$\endgroup$
  • $\begingroup$ This is a very good answer, but I would not like to call fvy[t] an exact solution since it is dependent on inexact parameters such as .0.001 and c. $\endgroup$ – m_goldberg Mar 24 '16 at 4:58
  • $\begingroup$ Yeah!Definately a good answer in this question ,but as a matter of fact I don't know the meaning of vy /. % /. C[1] -> cand FindRoot[%[0] == v0y, {c, .0016, 0, .00198}] (* {c -> 0.0019741} *) fvy = %% /. %may because of not expert in Mathematica,so excuse me can you explain a little nore clear?And also in my queston when ` k = 0.1` you give me the perfect fitting but when k=1,k=10.,ther will be some deviation,and also I was interested in this point.Thanks! $\endgroup$ – dcydhb Mar 24 '16 at 13:14
  • $\begingroup$ @m_goldberg You are correct. I have amended my overly generous description of fvy. Thanks. $\endgroup$ – bbgodfrey Mar 24 '16 at 18:53
  • $\begingroup$ @dcydhb I have added some explanatory material to my earlier answer. In addition, I added a new section to accommodate arbitrary constants. Perhaps, you will find cs = s[[2, 0, 1]] perplexing. I can only say that a combination of experience plus trial-and-error allowed my to choose the particular Part indices in this expression. There may well be a better way to obtain cs. $\endgroup$ – bbgodfrey Mar 24 '16 at 18:57
  • $\begingroup$ @bbgodfrey I have added several questions just now in my question and can you help me with that ,thanks a lot! $\endgroup$ – dcydhb Mar 26 '16 at 15:11
2
$\begingroup$

You "can't get the result" because data contains point {0, 282.843} while model contains a Log[t] term. Log[0] makes no sense.

To fix your code, you can shift the data at $t=0$ a little to e.g. $10^{-6}$:

data2 = {{10^-6, 282.843}, {1, 18.071}, {2,  6.2538}, {3, -3.58047}, 
        {4, -12.6429}, {5, -18.3376}};
fit2 = FindFit[data2, model, {a, b, c}, t];
Plot[Evaluate[model /. fit2], {t, 0, 5}, Epilog -> Point/@ data2]

enter image description here

Not bad at least for the first several points.

Also, you can remove the first point in data:

fit3 =FindFit[data[[2 ;;]], model, {a, b, c}, t]
Plot[Evaluate[model /. fit3], {t, 0, 30}, Epilog -> Point /@ data]

enter image description here

The fitting result isn't that good, but you've got a result this time. How to find a good fit is another story, and bbgodfrey has already given a good answer.

$\endgroup$
  • $\begingroup$ You can leave a comment right under the post you are referring to; Posting a comment as an answer is not the way to do it. $\endgroup$ – Sektor Mar 28 '16 at 12:13
  • $\begingroup$ ok.it's my first time write something in this community. I want to post a comment in the above answer initially ,but I don't have enough reputation.I have re-edited my answer. $\endgroup$ – jiaoeyushushu Mar 29 '16 at 3:28
  • $\begingroup$ @Sektor The original post is indeed an answer (of course a bit too brief). Now it should meet the standard. $\endgroup$ – xzczd Mar 29 '16 at 4:40

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.