Re: Help with solving ODE
- To: mathgroup at smc.vnet.net
- Subject: [mg81796] Re: Help with solving ODE
- From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
- Date: Wed, 3 Oct 2007 06:14:37 -0400 (EDT)
- Organization: The Open University, Milton Keynes, UK
- References: <fdvdq7$sul$1@smc.vnet.net>
Pioneer1 wrote: > Hi, > > Can anyone help solve this linearized differential equation: > > Iy'' + ky' = 2GMmd/a^2 > > Primes are time derivates of y (=theta=excursion angle). Is it > possible to solve this for the initial conditions y(0)=0 and y'(0)=0? > > I got the solution at sci.math for the non-linear version and I want > to compare the two. Here's the link to sci.math thread: > > http://groups.google.com/group/sci.math/browse_thread/thread/a6ee2f782df09625/53cf5573d354a3ab#53cf5573d354a3ab > > Further information is also available at sci.physics.research > > http://groups.google.com/group/sci.physics.research/browse_thread/thread/d391940cc173f9dc/eed90e6c3fee0edc#eed90e6c3fee0edc > > Parameters are: > >> y = theta = excursion angle in radians >> A = I = moment of inertia = 13,138,117.34 g cm^2 >> B = R = damping = for now I assume this to be zero >> C = k = torsion constant = 724.68 g cm^2 sec^-2 >> d = moment arm = 93.09 cm >> D = 2GMmd = 2 * 6.67*10^-8 * 158100 * 729.8 * 93.09 = 1432.82 >> a = distance between weights = 22.10 cm > > I would truly appreciate help with this. Thanks With Mathematica, you can use DSolve for an analytic solution. (In the example below, I have replaced capital 'I' by 'i' since capital I stands for the imaginary unit, Sqrt[-1], in Mathematica.) In[1]:= sol = DSolve[{i y''[t] + k y'[t] == 2 G M m d/a^2, y[0] == 0, y'[0] == 0}, y, t][[1]] Out[1]= (k t)/i (k t)/i 2 d G m M (i - E i + E k t) {y -> Function[{t}, -----------------------------------------]} (k t)/i 2 2 E (a k ) In[2]:= i y''[t] + k y'[t] == 2 G M m d/a^2 /. sol // Simplify Out[2]= True In[3]:= y[t] /. sol /. {i -> 13138117.34, k -> 724.68, G -> 6.67*10^(-8), M -> 158100, m -> 729.8, d -> 93.09, a -> 22.10} Out[3]= -6 7 7 0.0000551586 t (5.58621 10 (1.31381 10 - 1.31381 10 E + 0.0000551586 t 0.0000551586 t 724.68 E t)) / E In[4]:= % /. t -> 0 Out[4]= 0. In[5]:= D[%%, t] /. t -> 0 Out[5]= 0. In[6]:= $Version Out[6]= "6.0 for Microsoft Windows (32-bit) (June 19, 2007)" Regards, -- Jean-Marc