3 views

Skip to first unread message

Sep 15, 2007, 6:20:59 PM9/15/07

to

Hi,

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

Sep 15, 2007, 8:13:19 PM9/15/07

to

On 9 16 , 6 20 , Pioneer1 <1pione...@gmail.com> wrote:

> Hi,

>

> Can you please help obtain a numerical solution for this

>

> Ay'' + By' + Cy = D/(a - y)^2

>

> 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

Message has been deleted

Sep 15, 2007, 9:28:01 PM9/15/07

to

On Sep 15, 8:13 pm, bjoechan2...@gmail.com wrote:

>

> Is it then necessary to find the numerical solution?

>

> 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.

Sep 16, 2007, 3:49:38 AM9/16/07

to

Pioneer1 <1pio...@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.

--

Robert Israel isr...@math.MyUniversitysInitials.ca

Department of Mathematics http://www.math.ubc.ca/~israel

University of British Columbia Vancouver, BC, Canada

Sep 16, 2007, 10:35:19 AM9/16/07

to

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.

<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.

Sep 16, 2007, 5:28:43 PM9/16/07

to

Pioneer1 <1pio...@gmail.com> writes:

> 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.

Sep 18, 2007, 2:59:30 AM9/18/07

to

"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.

Message has been deleted

Sep 18, 2007, 8:44:41 AM9/18/07

to

On Sep 17, 2:28 am, Robert Israel

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

Message has been deleted

Sep 18, 2007, 11:05:00 PM9/18/07

to 1pio...@gmail.com

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.

<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.

Sep 19, 2007, 2:49:20 AM9/19/07

to

Pioneer1 <1pio...@gmail.com> writes:

> 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.

Sep 19, 2007, 9:37:16 PM9/19/07

to 1pioo...@gmail.com

On Sep 19, 2:49 am, Robert Israel

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.

Sep 20, 2007, 12:59:37 AM9/20/07

to

Pioneer1 <1pio...@gmail.com> writes:

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.

Sep 26, 2007, 8:57:13 PM9/26/07

to

On Sep 20, 12:59 am, Robert Israel

<isr...@math.MyUniversitysInitials.ca> wrote:

<isr...@math.MyUniversitysInitials.ca> wrote:

> 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

>

> 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?

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu