Non-linear ODE help
Pioneer1
Can you please help obtain a numerical solution for this
Ay'' + By' + Cy = D/(a - y)^2
The question arose in this thread
at sci.physics.research relating to the torsion pendulum motion.
y = theta = the excursion angle.
a = distance between weights
Usually the equation is linearized by eliminating y from the
denominator. I would like to find a solution with y. Can anyone help
obtain the numerical solution?
Thanks
bjoech...@gmail.com
> Hi,
>
> Can you please help obtain a numerical solution for this
>
> Ay'' + By' + Cy = D/(a - y)^2
>
If A, B, C, D and a are constants, I consider there is already an
analytical solution for it.
Is it then necessary to find the numerical solution?
best regards,
Chi Ho, Chan
Pioneer1
>
> Is it then necessary to find the numerical solution?
No. I was told that this did not have analytical solution. I would
appreciate it if you could let me know the analytical solution.
Many thanks.
Robert Israel
> Can you please help obtain a numerical solution for this
>
> Ay'' + By' + Cy = D/(a - y)^2
To obtain numerical solutions, you need to specify values for the
parameters A,B,C,D,a, an initial condition, and the values of t for
which you want the result.
--
Robert Israel isr...@math.MyUniversitysInitials.ca
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia Vancouver, BC, Canada
Pioneer1
<isr...@math.MyUniversitysInitials.ca> wrote:
> Pioneer1 <1pione...@gmail.com> writes:
> > Can you please help obtain a numerical solution for this
>
> > Ay'' + By' + Cy = D/(a - y)^2
>
> To obtain numerical solutions, you need to specify values for the
> parameters A,B,C,D,a, an initial condition, and the values of t for
> which you want the result.
Thanks. I list the constants below. But I made a change in the
equation, in the force equation y should be y*d:
Ay'' + By' + Cy = D/(a - y*d)^2
And I do not know the damping constant B yet, so the equation reduces
to
Ay'' + Cy = D/(a - y*d)^2
Constants are:
y = theta = excursion angle in radians
A = I = moment of inertia = 13,138,117.34 g cm^2
B = R = [not known at this point]
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
Initial condition:
y = 0 and t = 0
Range of t:
t = 0 to t = 840 seconds
The pendulum had a period of 420 seconds. But this was what we would
today call half period.
Thanks again for your help.
Robert Israel
> On Sep 16, 3:49 am, Robert Israel
> <isr...@math.MyUniversitysInitials.ca> wrote:
> > Pioneer1 <1pione...@gmail.com> writes:
> > > Can you please help obtain a numerical solution for this
> >
> > > Ay'' + By' + Cy = D/(a - y)^2
> >
> > To obtain numerical solutions, you need to specify values for the
> > parameters A,B,C,D,a, an initial condition, and the values of t for
> > which you want the result.
>
> Thanks. I list the constants below. But I made a change in the
> equation, in the force equation y should be y*d:
>
> Ay'' + By' + Cy = D/(a - y*d)^2
>
> And I do not know the damping constant B yet, so the equation reduces
> to
>
> Ay'' + Cy = D/(a - y*d)^2
It only reduces to that if B is 0. Unknown is not the same as 0.
> Constants are:
>
> y = theta = excursion angle in radians
> A = I = moment of inertia = 13,138,117.34 g cm^2
> B = R = [not known at this point]
> 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
>
> Initial condition:
>
> y = 0 and t = 0
For a second-order differential equation, you also need to specify y'(0).
> Range of t:
>
> t = 0 to t = 840 seconds
With B = 0 and the initial conditions y(0) = 0 and y'(0) = 0, Maple tells me
that y(840) = .518143785620031448e-4, y'(840) = -.481752705908163793e-5.
For a plot of y(t) over that range, see
<http://www.math.ubc.ca/~israel/pioneer.gif>
> The pendulum had a period of 420 seconds. But this was what we would
> today call half period.
Yes, it seems pretty close to that.
Michael Jørgensen
"Pioneer1" <1pio...@gmail.com> wrote in message
news:1189953319.9...@r29g2000hsg.googlegroups.com...
> On Sep 16, 3:49 am, Robert Israel
> <isr...@math.MyUniversitysInitials.ca> wrote:
>> Pioneer1 <1pione...@gmail.com> writes:
>> > Can you please help obtain a numerical solution for this
>>
>> > Ay'' + By' + Cy = D/(a - y)^2
>>
>> To obtain numerical solutions, you need to specify values for the
>> parameters A,B,C,D,a, an initial condition, and the values of t for
>> which you want the result.
>
> Thanks. I list the constants below. But I made a change in the
> equation, in the force equation y should be y*d:
>
> Ay'' + By' + Cy = D/(a - y*d)^2
>
> And I do not know the damping constant B yet, so the equation reduces
> to
>
> Ay'' + Cy = D/(a - y*d)^2
This equation is solvable using elliptic functions: Multiply by y' and
integrate on both sides with respect to t. This gives an equation of the
form:
(y')^2 = P3(y) / (a-y*d), where P3() is a cubic polynomial. The coefficients
of this polynomial are derived from the constants of the original equation,
as well as the initial conditions.
The solutions are periodic, and the period is expressed in terms of complete
elliptic integrals. The roots of the cubic polynomial determine the
inflection points, i.e. minimum and maximum value of y(t). Writing down the
analytical solution involves determining the roots of the cubic polynomial.
This becomes quite messy, unless you already know some of the roots.
If you e.g. know that the minimum value of y(t) is y=0, then that is a root
of the cubic. You may then write P3(y) = k1 * y * (y - y1) * (y - y2), where
k1, y1, and y2 are constants. Now a table of integrals will give you the
complete analytical solution in terms of these constants and a and d.
If the damping is very small, then you may be able to get some quantitative
results using perturbation theory. For instance, the undamped case has a
conserved quantity corresponding to the total energy of the system. You may
be able to estimate the rate of decrease of that energy over time, as caused
by the damping.
Hope that helps.
-Michael.
Narasimham
The non-linear dynamical oscillation should be quite interesting
from maths and physics viewpoints; to me it appears as a different
sort of
non-linear pendulum, where the forcing function is now gravitational.
At y = theta = a/d, infinite accelerations might force some different
mode of path retracing as the forcing function (torque)is always
positive. In your plot, y maximum value is far below this critical
angle ycrit = a/d =0.2376 radians.So I request you also upload plots
of
y-y' and y-y'' ovals for a few assumed damping numerical values of B
to see how gravity/damping bring about non-linearity.Right now,its
behavior is not so far away from a linear harmonic situation.
Regards,
Narasimham
Pioneer1
<isr...@math.MyUniversitysInitials.ca> wrote:
> With B = 0 and the initial conditions y(0) = 0 and y'(0) = 0, Maple tells me
> that y(840) = .518143785620031448e-4, y'(840) = -.481752705908163793e-5.
>
the damping constant but I assume that for one oscillation damping
doesn't have an effect.
But I am confused about the numbers. You have y(840) = .518*10^-4
(radians?)
This is the angle at t=840 seconds. I assume that this is same as
t=420 seconds
y(420) = .518*10^-4 radians
Is this correct?
For this experiment (experiment IV in Cavendish's original paper) I
calculated the maximum excursion as follows:
The arm moves from 31.3 divisions to 17.1 divisions. The "point of
rest" is given as 24.02 divisions.
>From this I compute that the maximum excursion is
31.3 - 24.02 = 7.28 divisions
Each division is 1/20 inches and 1 inch is 2.54 cm so,
7.28 * 1/20 * 2.54 = 0.9245 cm
Moment arm is 93.09 cm so,
0.9245/93.09 = 0.0099 radians
You have y(840) = 0.0000518 radians. This is about 200 times smaller
than Cavendish's value.
What are your thoughts on this? Why such a big discrepency? Many
thanks again for helping.
Robert Israel
> On Sep 16, 5:28 pm, Robert Israel
> <isr...@math.MyUniversitysInitials.ca> wrote:
> > With B = 0 and the initial conditions y(0) = 0 and y'(0) = 0, Maple tells
> > me
> > that y(840) = .518143785620031448e-4, y'(840) = -.481752705908163793e-5.
> >
> Thank you very much for your help. This is great! I will try to find
> the damping constant but I assume that for one oscillation damping
> doesn't have an effect.
>
> But I am confused about the numbers. You have y(840) = .518*10^-4
> (radians?)
If your original ODE had y in radians, yes.
> This is the angle at t=840 seconds. I assume that this is same as
> t=420 seconds
No. Didn't you look at the plot? And why would you assume that?
The period is approximately 840, there are minima at t=0 and approximately
840, and a maximum at approximately t=420.
> y(420) = .518*10^-4 radians
> Is this correct?
No, y(420) looks to be a bit more than .008.
Pioneer1
Ok. Thanks. I looked at the plot I don't know how I didn't see this.
0.008 looks closer to what I calculated. Would it be possible to get
the data so that I can look at it in excel and maybe compare it with
the linear solution. The plot looks linear but taking residuals may
give a better idea. Thanks again.
Robert Israel
OK, here are the values of t, y(t) and y'(t) for t from 0 to 870 by 10:
0, 0.00000000, 0.00000000
10, 0.00001116, 0.00000223
20, 0.00004458, 0.00000445
30, 0.00010006, 0.00000665
40, 0.00017730, 0.00000881
50, 0.00027590, 0.00001092
60, 0.00039535, 0.00001297
70, 0.00053501, 0.00001496
80, 0.00069417, 0.00001687
90, 0.00087197, 0.00001868
100, 0.00106747, 0.00002040
110, 0.00127961, 0.00002201
120, 0.00150725, 0.00002350
130, 0.00174918, 0.00002487
140, 0.00200412, 0.00002610
150, 0.00227070, 0.00002720
160, 0.00254753, 0.00002815
170, 0.00283312, 0.00002895
180, 0.00312595, 0.00002959
190, 0.00342447, 0.00003008
200, 0.00372707, 0.00003041
210, 0.00403216, 0.00003058
220, 0.00433813, 0.00003059
230, 0.00464334, 0.00003043
240, 0.00494617, 0.00003011
250, 0.00524500, 0.00002963
260, 0.00553825, 0.00002899
270, 0.00582436, 0.00002820
280, 0.00610183, 0.00002726
290, 0.00636917, 0.00002618
300, 0.00662495, 0.00002496
310, 0.00686783, 0.00002360
320, 0.00709651, 0.00002212
330, 0.00730979, 0.00002052
340, 0.00750652, 0.00001881
350, 0.00768567, 0.00001700
360, 0.00784627, 0.00001510
370, 0.00798749, 0.00001313
380, 0.00810858, 0.00001108
390, 0.00820889, 0.00000897
400, 0.00828788, 0.00000682
410, 0.00834514, 0.00000463
420, 0.00838036, 0.00000241
430, 0.00839336, 0.00000019
440, 0.00838408, -0.00000204
450, 0.00835256, -0.00000426
460, 0.00829896, -0.00000646
470, 0.00822356, -0.00000862
480, 0.00812678, -0.00001073
490, 0.00800912, -0.00001279
500, 0.00787121, -0.00001478
510, 0.00771379, -0.00001669
520, 0.00753769, -0.00001852
530, 0.00734383, -0.00002024
540, 0.00713324, -0.00002186
550, 0.00690704, -0.00002336
560, 0.00666643, -0.00002474
570, 0.00641270, -0.00002599
580, 0.00614720, -0.00002709
590, 0.00587132, -0.00002806
600, 0.00558654, -0.00002887
610, 0.00529435, -0.00002954
620, 0.00499632, -0.00003004
630, 0.00469405, -0.00003039
640, 0.00438912, -0.00003057
650, 0.00408317, -0.00003059
660, 0.00377780, -0.00003045
670, 0.00347465, -0.00003015
680, 0.00317532, -0.00002969
690, 0.00288143, -0.00002907
700, 0.00259452, -0.00002829
710, 0.00231612, -0.00002736
720, 0.00204770, -0.00002629
730, 0.00179071, -0.00002508
740, 0.00154650, -0.00002374
750, 0.00131639, -0.00002227
760, 0.00110158, -0.00002068
770, 0.00090323, -0.00001898
780, 0.00072238, -0.00001718
790, 0.00056001, -0.00001528
800, 0.00041698, -0.00001331
810, 0.00029405, -0.00001127
820, 0.00019186, -0.00000916
830, 0.00011097, -0.00000701
840, 0.00005181, -0.00000482
850, 0.00001470, -0.00000260
860, -0.00000016, -0.00000037
870, 0.00000728, 0.00000186
As you can see, the period actually seems to be about 860, not 840.
Pioneer1
<isr...@math.MyUniversitysInitials.ca> wrote:
>
> 0, 0.00000000, 0.00000000
Many thanks for this. I entered the values to excel but I realized
that the linear solution too must be solved for the same initial
conditions. Can you help with that as well so that I can compare them?
The linearized differential equation is:
Iy'' + ky' = 2GMmd/a^2
Primes are time derivates of theta as before. Is it possible to solve
this for the initial conditions y(0)=0 and y'(0)=0?