Message Boards Message Boards

0
|
4747 Views
|
10 Replies
|
0 Total Likes
View groups...
Share
Share this post:

Numerical solving of second order differential equation?

Marc van Gestel
Posted 11 years ago
Hi all,

New to the community so I might be listing in the wrong groups. 
I'm trying to fathom the mathematics of the working of a Trebuchet, a giant medieval slingshot.

The problem is this: I have a second derivative Y'' = - C * Sin(Y)
Y is dependent on t (time). y(0) = 3/4 * Pi and y'(0) = 0. What I'd like to know is how to phrase the whole thing in WolfranAlpha if I want to solve what t is at y(t) = Pi/4?

C = 26.48
I arrived at: {y'' = - 26,48 * Sin (y), y'(0) = 3/4*Pi, y'(0) = 0}
But hten I only get the solution y and no other thing but a plot. I need to interact with values.....

How do I state that I want to know t at y = Pi/4?

thanks for your replies.....

Marc
POSTED BY: Marc van Gestel
10 Replies
Sort By:
Markus Schopfer
Markus Schopfer
Posted 11 years ago
Not sure if I understood your problem correct, but maybe tha's what you could use:
c = 26.48;
val = Pi/4;

FindRoot[(y[t] /. 
    DSolve[{y[0] == 3/4*Pi, y'[0] == 0, y''[t] == -c*Sin[t]}, y,
     t]) == val,
{t, 1}]
 
POSTED BY: Markus Schopfer
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Hi and thanks. You understand correctly but if I try to copy-paste this into the window at WWW.WolframAlpha.com it doesn't recognize the syntax. As a program in some compiler it would work I think.
Marc
POSTED BY: Marc van Gestel
Markus Schopfer
Markus Schopfer
Posted 11 years ago
Ah, i overlooked that you refer to Wolfram Alpha.
My solution is written in Mathematica-code.

I don't know if you can express this in one phrase foe wolfram alpha.
However you can of course you can do it in two steps.
Just first solve the differential equations:
Differential equation solution {y(0)=3/4*Pi,y'(0)=0,y''(t)=-26.48*sin(t)}}}

then press the solution:
y(t) = -26.48 t+26.48 sin(t)+2.35619

add "=Pi/4" and execute this query again.

If you use Wolfram Aplha pro, you can of course add "=x" instead and enable interactivity.
 

 
POSTED BY: Markus Schopfer
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Hi, your suggestion works great, Thanks!
Being curious now as to what Mathematica is and how it works (I used to code a lot in Delphi and Basic) I donloaded a trial version for a month. But how do I start with your code in Mathematica?
(excuse me if I bother you with these trivial questions....)
thnx, Marc
POSTED BY: Marc van Gestel
Markus Schopfer
Markus Schopfer
Posted 11 years ago
Hi Marc, you're not bothering me.
The community sites are here for asking questions.


As a programmer, you will love Mathematica.

To execute my code, simply open a mathematica notebook and paste the code into it.
Press Shift+Enter and you will get the answer (in the form of a rule): {t -> 0.714745}.

Or try this one:

 (* Just in case we have stored something in these variables, we clear \
 them *)
 Clear[f, c, y, t, x];
 
 (* Let Mathematica Solve the Differetial equation y(t).
 In case we have more solutions, we use only the first one (with First).
 Since we know y(t) and want to get t, we are looking for the inverse \
 function of y(t) *)
 
f[c_, x_] =
  InverseFunction[(y /.
      First[DSolve[{y[0] == 3/4*Pi, y'[0] == 0, y''[t] == -c*Sin[t]},
        y, t]])][x] ;


(* Finally we plot t=f[c,x] *)
Manipulate[
Plot[f[c, x] , {x, -Pi, 2 Pi}, PlotRange -> {{-Pi, 2*Pi}, {-2, 2}},
  AxesLabel -> {"x", "t"}],
{{c, 26.48}, -50, 50, Appearance -> "Open"}]

Just play around a bit with mathematica and work through the documentation and tutorials.
It's an incredibly powerful tool. 
POSTED BY: Markus Schopfer
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Hi Markus,

This is working like a charm! Thanks a lot.....working throught the manuals right now.

Marc
POSTED BY: Marc van Gestel
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Hi Markus,

Digging in the documentation for Mathematica. Just to make sure.... how would I rephrase your notebookfile if I wanted to search the function y(t)'' = -26.48*sin(y(t))? This is the original function.... Your notebookfile regards the function y'' - -26.48*sin(t)

regards,

Marc
POSTED BY: Marc van Gestel
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Hi Markus,

I found the next notebookfile but it is not displaying the Plots. How come?
thx for your time....

Marc
 (* Solve the seesaw treb problem *)
 (*seesaw1.0.nb DB Siano Aug. 21, 1997*)
 m1=100;m2=1;l1=1;l2=4;g=32;(* the parameters*)
 
 c=g*(l1 m1 - l2 m2)/(m1 l1 l1+ m2 l2 l2);
 
 (*Solve the DE:*)
 solss=NDSolve[{th''[t]==-c Sin[th[t]],th[0]==3 Pi/4,th'[0]==0},{th[t]},{t,0, 1}];
 (* get th and the velocity from this solution*)
thint[t_]=Chop[th[t]/.Flatten[solss][[1]]];
v[t_]=l2*thint'[t];(*tangential velocity of the projectile*)

Print["time when th=pi/4 is ",tsolss=t/.FindRoot[thint[t]==Pi/4,{t,.2,.4}]];
Print["vel at th=pi/4 is=",v0pi4=l2*thint'[tsolss]];
Print["range when release is at Pi/4 is ",r=v0pi4^2/g," ft"];
Print["thoret max range=",rth=N[2 m1 l1 (1+Sin[3 Pi/4])/m2]," ft"];
Print["efficiency=",r/rth];


(* Plot the results*)
Plot[{thint[t]},{t,0,.5},PlotLabel->"th[t]"];
(* plot the range as a function of the time at release*)
range[t_]=2*v[t]*v[t]*Cos[thint[t]]*Sin[thint[t]]/g;
Plot[range[t],{t,0,.5},PlotLabel->"range[t]"];
POSTED BY: Marc van Gestel
Markus Schopfer
Markus Schopfer
Posted 11 years ago
Hi Marc.
Just remove the semicolons (;) from the end of the Lines with the Plot[..] ; commands.
  (* Solve the seesaw treb problem *)
  (*seesaw1.0.nb DB Siano Aug. 21, 1997*)
  m1=100;m2=1;l1=1;l2=4;g=32;(* the parameters*)
 
  c=g*(l1 m1 - l2 m2)/(m1 l1 l1+ m2 l2 l2);
 
  (*Solve the DE:*)
  solss=NDSolve[{th''[t]==-c Sin[th[t]],th[0]==3 Pi/4,th'[0]==0},{th[t]},{t,0, 1}];
  (* get th and the velocity from this solution*)
thint[t_]=Chop[th[t]/.Flatten[solss][[1]]];
v[t_]=l2*thint'[t];(*tangential velocity of the projectile*)

Print["time when th=pi/4 is ",tsolss=t/.FindRoot[thint[t]==Pi/4,{t,.2,.4}]];
Print["vel at th=pi/4 is=",v0pi4=l2*thint'[tsolss]];
Print["range when release is at Pi/4 is ",r=v0pi4^2/g," ft"];
Print["thoret max range=",rth=N[2 m1 l1 (1+Sin[3 Pi/4])/m2]," ft"];
Print["efficiency=",r/rth];


(* Plot the results*)
Plot[{thint[t]},{t,0,.5},PlotLabel->"th[t]"]
(* plot the range as a function of the time at release*)
range[t_]=2*v[t]*v[t]*Cos[thint[t]]*Sin[thint[t]]/g;
Plot[range[t],{t,0,.5},PlotLabel->"range[t]"]

A semicolon alway supresses the output. 
POSTED BY: Markus Schopfer
Marc van Gestel
Marc van Gestel
Posted 11 years ago
Great! Simple as that......
POSTED BY: Marc van Gestel
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract