Help solving linear ODE
Pioneer1
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:
Further information is also available at sci.physics.research
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
Jean-Marc Gulliet
Here is the analytic solution with Maple.
|\^/| Maple 11 (IBM INTEL NT)
._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc.
2007
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> dsolve([I*((D@@2)(y))(t)+k*(D(y))(t) =\
> 2*G*M*m*d/a^2, y(0) = 0, (D(y))(0) = 0]);
bytes used=4000568, alloc=3276200, time=0.50
2 I exp(k t I) G M m d 2 G M m d t 2 I G M m d
y(t) = ---------------------- + ----------- - -----------
2 2 2 2 2
k a a k k a
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
Pioneer1
wrote:
> Pioneer1 wrote:
>
> 2 I exp(k t I) G M m d 2 G M m d t 2 I G M m d
> y(t) = ---------------------- + ----------- - -----------
> 2 2 2 2 2
> k a a k k a
>
Thank you very much for your help with this. This was very helpful.
I have one question. When I entered the formula in excel exp(ktI)
gives #NUM error message. So I used the Mathematica version which has
exp(kt/I). Do you know why the discrepency?
Thanks again for helping.
Jean-Marc Gulliet
Hum, yes (shame on me). Indeed, it's my mistake: I forget to change the
uppercase I into a lowercase one in Maple (and obviously did not check
the result). In Maple too, capital I stands for the imaginary unit
sqrt(-1). So, this time with the correct syntax, Maple 11.01 returns
> dsolve([i*((D@@2)(y))(t)+k*(D(y))(t) = 2*G*M*m*d/a^2, y(0) = 0,
(D(y))(0) = 0]);
/ k t\
2 i exp|- ---| G M m d
\ i / 2 G M m d t 2 G M m d i
y(t) = ---------------------- + ----------- - -----------
2 2 2 2 2
k a a k k a
Also, tweaking Mathematica a little bit, we can check that the solutions
returned by both systems are equivalent.
In[1]:= y[t] /. DSolve[{i y''[t] + k y'[t] == 2 G M m d/a^2, y[0] == 0,
y'[0] == 0},
y, t][[1]] // ComplexExpand
Out[1]=
2 d G i m M 2 d G i m M 2 d G m M t
-(-----------) + ---------------- + -----------
2 2 (k t)/i 2 2 2
a k E (a k ) a k
Sorry about the confusion.
Regards,
--
Jean-Marc
Pioneer1
wrote:
>
> / k t\
> 2 i exp|- ---| G M m d
> \ i / 2 G M m d t 2 G M m d i
> y(t) = ---------------------- + ----------- - -----------
> 2 2 2 2 2
> k a a k k a
>
thought that if I can have y'' and y' and substitute those in the
differential equation I should be able to solve for R. Can you obtain
some value of y' and y'' so that I can solve for R?
The differential equation is:
iy'' + Ry' + ky = 2GMmd/a^2
Thanks again for your help.
Jean-Marc Gulliet
Here is what is returned by Mathematica 6.0.1. (Sorry, I am not near a
machine with Maple 11.01 installed on it. I must wait until tomorrow to
access it.)
(* We built a list of replacement rules for the known parameters. I have
kept the values given in your first post and also the initial conditions. *)
In[1]:= params = {i -> 13138117.34, k -> 724.68, G -> 6.67*10^(-8),
M -> 158100, m -> 729.8, d -> 93.09, a -> 22.10};
(* We solve the ODE *)
In[2]:= sol =
DSolve[{i y''[t] + R y'[t] + k y[t] == 2 G M m d/a^2, y[0] == 0,
y'[0] == 0}, y, t][[1]]
(* Note that the ASCII string "\[ExponentialE]" (without the double
quotes) is a symbol for the nice display of the exponential function E. *)
Out[2]= {y ->
Function[{t}, -(1/(a^2 k Sqrt[-4 i k + R^2]))
d G m M (-\[ExponentialE]^((((-R - Sqrt[-4 i k + R^2]) t)/(2 i)))
R + \[ExponentialE]^(((-R + Sqrt[-4 i k + R^2]) t)/(2 i))
R - 2 Sqrt[-4 i k +
R^2] + \[ExponentialE]^(((-R - Sqrt[-4 i k + R^2]) t)/(2 i))
Sqrt[-4 i k +
R^2] + \[ExponentialE]^(((-R + Sqrt[-4 i k + R^2]) t)/(2 i))
Sqrt[-4 i k + R^2])]}
In[3]:= i y''[t] + R y'[t] + k y[t] == 2 G M m d/a^2 /.
sol // Simplify
Out[3]= True
(* Here is the first derivative *)
In[4]:= D[y[t] /. sol, t] // FullSimplify
Out[4]= (2 d (-1 + \[ExponentialE]^((Sqrt[-4 i k + R^2] t)/
i)) G m M)/(\[ExponentialE]^(((R + Sqrt[-4 i k + R^2]) t)/(
2 i)) (a^2 Sqrt[-4 i k + R^2]))
In[5]:= % /. params
Out[5]= (2.9336617706522796 (-1 + \[ExponentialE]^(
7.61144*10^-8 Sqrt[-3.80837*10^10 + R^2] t)))/(\[ExponentialE]^(
3.80572*10^-8 (R +
Sqrt[-3.80837*10^10 + R^2]) t) Sqrt[-3.80837*10^10 + R^2])
(* Here is the second derivative *)
In[6]:= D[y[t] /. sol, {t, 2}] // FullSimplify
Out[6]= (d G m M (R - \[ExponentialE]^((Sqrt[-4 i k + R^2] t)/i)
R + (1 + \[ExponentialE]^((Sqrt[-4 i k + R^2] t)/
i)) Sqrt[-4 i k + R^2]))/(\[ExponentialE]^(((R +
Sqrt[-4 i k + R^2]) t)/(2 i)) (a^2 i Sqrt[-4 i k + R^2]))
In[7]:= % /. params
Out[7]= (1.11647*10^-7 (R - \[ExponentialE]^(
7.61144*10^-8 Sqrt[-3.80837*10^10 + R^2] t)
R + (1 + \[ExponentialE]^(
7.61144*10^-8 Sqrt[-3.80837*10^10 + R^2]
t)) Sqrt[-3.80837*10^10 + R^2]))/(\[ExponentialE]^(
3.80572*10^-8 (R + Sqrt[-3.80837*10^10 + R^2]) t)
Sqrt[-3.80837*10^10 + R^2])
Regards,
--
Jean-Marc