[Date Index] [Thread Index] [Author Index]

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