Beschrijving van de rotaties van een star lichaam met Eulerhoeken of Eulerparameters?

(1)

Beschrijving van de rotaties van een star lichaam met

Eulerhoeken of Eulerparameters?

Citation for published version (APA):

Deen, P. (1986). Beschrijving van de rotaties van een star lichaam met Eulerhoeken of Eulerparameters? (DCT rapporten; Vol. 1986.056). Technische Hogeschool Eindhoven.

Document status and date: Gepubliceerd: 01/01/1986

Document Version:

Uitgevers PDF, ook bekend als Version of Record

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

(2)

(3)

2

7 '7

'7

1

9

2

- i

-

(4)

(5)

- -.,,I ... ... _.___I .._ ___. _ _ . . . . .. _.-.- .... .

-

. . .

_-

. . . . . . . . .

__

... _ _ _ _ _ ~ I.-_. .... . ... . . , . . . . . .

__

...

__

. _l___l__ ... -. . -L.- . . . . . . . . .... > ... " . . . .I ._I.I . . . . _I . . . . . ~~ .... . . . * . . . .~ - . __ . . . I ... ._I___I.____.____^.-- ... .... .. . . . . .- "

-3

-

(6)

(7)

(8)

(9)

(10)

(11)

( V n =

.u PLI

h

Y ' N

(12)

de

(13)

(14)

(15)

(16)

(17)

'Z, r.2 a i 10 I

-

_

-

_

-2

o

0

Q. I

-

-0.3

3 - - - -

ì u

c:

o

_ -

- -i .

c

iw

_ _ - - -

i o

I

_ -

e

(18)

(19)

(20)

(21)

(22)

I

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

c

L"

c

(32)

(33)

L

c

e

c

e?

c

c? L

G

P ec

c

(34)

(35)

i f

c

<

c

<

{

a

c

6

<

6

<

e

t f k

c

d

c

i I

a

c

e

f

c

e

a

e

c

I

e

d

f t

<

d

$

(36)

(37)

(38)

c

G et

e

c

G

c

>- e_,

(39)

(40)

(41)

(42)

(43)

(44)

(45)

(46)

(47)

óQ

16

T P 1

(48)

(49)

(50)

i.

c

(51)

(52)

<

i f

c

f

c

<

c

d

c

f i' b'

c

E E { b' i i

c

i l

c

i i: i i

a

c

i: E

<

c

i

c

E

a

i

<

f i i' i 4

e

c

5 t

(53)

(54)

(55)

(56)

(57)

(58)

(59)

(60)

EULERHOEKEN . .

I

-

EULERPARAMETERS ,001 I 1 .OO 1.49 1.99 2.49 2.99 .OO ,50 TIJD ( = X )

-

. _ _ _ _ _ _ _ ._ - - . -. HKNf01 :1 .100E 2 :2 .l00E 2 I3 ,100E 2 11 ,100E 1 12 .200E 1 13 .300E 1 -OL ,100E -8 :PU-TIME .G69E PC I THETA FI HKN/01 I1 ,100E 2 I2 .100E 2 .100E 2 I3 I 1

.ImE

1 .200E 1 12 13 .300E 1 TOL .ímE -8 :PU-TIME .G69E 1 Q0 P1 EULERHOEKEN /

/

/’.

/ EULERPARAMETERS /’ ,/‘ ,/“ TIJD ( = X ) HKN/0í [ 1 .100E 2 [2 .í00E 2 [3 ,100E 2 11 .100E 1 12 .200E 13 I 300E ‘OL .í00E :PU-TIME .G69E OMEGA1 OMEGA2 ._-.-_ OMEGA3 11 .100E 2 12 .100E 2 13 .100E 2 11 .100E 1 12 .200E 1 13 .3@0E 1 roL ,100E -8 :PU-TIME .G69E 82 1 83

(61)

I

(62)

1

PAR/01 EULERPARAMETERC .100E 2 ,100E 2 .100E 2 .Í00E 1 M2 .200E 1 M3 .300E 1 TOL .100E -8 CPU-TIME .402E 1 QO EULERPARAMETERC TIJD (=X) EULERPARAMETERS /

/

/’ / ,,“ ,_,- /// ,,/’ /’ ,/-

/

,,.f. ,/-- ,/‘ ,/’ ,/

/

/’ ,,/- ,/’

/

/’

, 1 .50 1 .OO 1.50 2.00 2.50 3.00 TIJD (=X) PAR/01 11 .100E 2 12 ,100E 2 13 ,100E 2 Y1 .100E 1 Y2 .200E 1 Y3 .300E 1 TOL .Í00E -8 -PU-TIME .402E Q2 1 PAW01 ! I ,100E 2 12 .100E 2 [3 ,100E 2 11 .100E 1 12 .200E I 13 .300E I D L .100E -8 :PU-TIME ,402E 1 OMEGA1 OMEGA2 OMEGA3

(63)

b m < PARAMETERS PARAMETERS & -5 PARAMETERS

(64)

7.n

FARAMETERS

PARAMETERS

W

(65)

I N39N ISlVVld83h NNNmmmmN N39NISlVVld83A

49.

22 j 8/ N t3 i- &&.?

(66)

HKNl02DET EULERHOEKEN

l

.83.

'

O

Q

~

2 .a3 1 .@0 0-

'

:J.cz2 26 HKN/02DET I1 .100E 2 I2 .100E 2 I3 .180E 2 M1 J00E 8 M2 .%BE 0 M3 .!300E 0 TOL

.lmE

-a

CPU-TIME .1%E 2 IiKN/02DEí EULERHOEKEN 11 .188E 2 12 .180E 2 13 .100E 2 MI .000E Q M2 .800E 0 M3 .80!3E 0 TOL .100E -8 CPU-TIME .í93E 2 OMEGA1 OMEGA2 .___.___.._..__._._....- I

.

a01 1 ~

.OO .i7 .33 .5a .67 -83 1.80

OMEGA3 TIJD (=X) X1Q-'rz.i222C

~ 3 U Z Z Z C

(67)

!C

..!?!a-

& I P

w

< 3 x z 6) i kl m PARAMETERS PARAMETERS

(68)

-_

EULERPAR AMETERC /' .-.a0. .30 .83 í .G7 2 5 0 TIJD - 1 3.33 4.17 5 3 0 X ) PAR/02 1 .í0OE 2 2 ,100E 2 3 .10PE 2 Ií .000E 0 I2 .000E 0 i3 .B0QE 0 'OL ,100E -8 :PU-TIME .580E í QO Q1 ---- ~ , PAR/02 TIJD -.33.. -.M.. /,' /,"

/"

//' //' i /" /. i 3.33 4.17 5-08 =X ) FULERPARAMETERS í .00, 1 .a0 .83 í .F7 2.50 3.33 4.17 5.00 TIJD i-X) :í .í00E 2 :2 .í00E 2 : 3 ,100E 2 11 .0@@E 0 12 J00E 0 13 300E 0 TOL ,100E -8 :PU-TIME .580E Q2 Q3 1 PARI02 !í .î0OE 2 12 .í0@E 2 13 ,100E 2 11 .000E 0 12 .@@@E 0 13 .000E 0 raL . í m ~ -8 :PU-TIME .5WE í OMEGAI OMEGA2 OMEGA3

(69)

+- N w 4 lx In i ? N ?

c

e I m m +- N m U e

a:

? I m W 8 \ ?

(70)

W m h

.-.

U v Zi -/ u; -I NI

T

n

-

U Y

z

SSi313WVdVd AI

zj

8

T

(71)

.r

"\ x101 EULERHOEKEN . l O l .21- ,10- ____-_-- _____ _____--- - --- ______._.-- -_- . .._.___ --- I .50 2.50 2.50 3,.00 TIJD (=X) EULERPAR AMETERS

'

0

1

I 3.00 1 .OO 2.50 HKN/IIII !I .îQ0E 2 !2 .î00E 2 [3 .100E 2 I 1 -.100E 1 12 ,000E 0 13 .000€ 0 roL .1 OOE -8 :PU-TIM€ .190E 2 PSI THETA FI HKN/03 1 .1@0E 2 2 .100€ 2 3

.

l00E 2 I1 -.I#E 1 12 ,000E 0 13 .000E 0 'OL .100€ -8 :PU-TIME ,190E 2 00 Q1 EULERHOEKEN XI 0-1 -.5Q 1 .50 1-00 1.50 2.00 2.50 3.0Q EULERHOEKEN XI 0-1 -.5Q. 1 2.00 2.50 3.0Q TIJD (=X) Lu u> I- 4 4 -1 b %I .50- I-

L.50~

\\\\\\

-2.00- -2.50.. -3.00- HKN/03 1 .1QQE 2 2 .í00E 2 3 .100E 2 1 -.100E 1 2 .000E 0 3 .000E 0 OL .180E -8 PU-TIN€ .190E 2 OMEGA1 OMEGA2 OMEGA3 HKN/03 [ I .100E 2 12 .I00E 2 13 .í00E 2 11 -.100E 1 12 ,000E 0 13 .000E 0 roL

.

I OOE -8 :PU-TIME .190E 2 Q2 Q3

(72)

EULERPARAMETERC

-

'221

.-_._ , 1 1.80 i:5B---- 2 .88 2.50 3.80 PAR/03 ' 1 .í00E 2 '2 .100E 2 '3 .188E 2 11 -,l@@E 1 12 .880E 0 13 .080E 0 .OL .I00E -8 :PU-TIME .212E 1 P0 PI EULERPARAMETERC x 1 0-5 I FULERPARAMETERC xlO-1 1 1.50 2.80 2.50 3.88 TIJD ( = X ) # E5 c.-1

.m-

2

d

a

PAR103 11 .lQ0E 2 12 .10QE 2 13 J00E 2 M1 -,IEfJE 1 Y2 .80QE 0 13 .800E 8 TOL .100E -8 :PU-TIME .212E 1 P2 PAR/03 11 .100E 2 12

.

l00E 2 13 .100E 2 M1 -.íQ8E 1 M2 J00E @ M3 .880E 0 TOL .í00E -8 IPU-TIME .212E 1 OMEGA1 OMEGA2 .- .--.. _.._..._ ^ _ OMEGA3

(73)

m m 8

b-

8

_W o. z Y I I m m n AI U o

-

Nj ij

vi

m m m m

b-

e W m Q1 I Y 1 -$LJ 5-

2

(74)

d

HKN/04 I .10@E 2 2 ,100E 2 3

.

l08E 2

I1

,100E 2 12 ,000E 0 13 .800E @ ’OL . I 00E -8 :PU-TIME .752E 2 PSI THETA FI HKN44 FULERPARAMETERS I1 ,100E 2 I2 . I ~ E 2 I3 .I00E 2 MI

.ímE

2 M2 .000E 0 M3 .B@@E 0 TOL .10@E -8 CPU-TIME .752E 2 nE, &- GI Is I I EULERHOEKEN 4.99- .m 1

.

00 I83 1.66 2.50 3.33 4.16 4.99 TIJD ( = X ) EULERPARAMETERC x1 0-5 .G7_ v) a 1

.

,0 .83 I .GG 3 . 3 TIJD (=X) -.33. -.G& IX##

V

HKN/04 I í .100E 2 12 .100E 2 13 .í00E 2 Y 1 .ImE 2 Y2 .@BOE 0 Y3 .0@@E 0 TOL ,108E -8 IPU-TIME .I52E 2 OMEGAI OMEGA2 OMEGA3 HKNl04 11 .í00E 2 12 .l@@E 2 13 .í0@E 2 Y1 ,í00E 2 Y2 .000E 0 Y3 .088E 0 TOL ,100E -8 >PU-TIME .752E 2 Q3

(75)

in m Lu I-.

z

8 Gj

(76)

$9.

li o) m Sti313liVP1Vd

i

ai.

62

(77)

O02

-

Ordinary Differemkl bqu:;tfons

W2BDF

NOTE before wing this rctl

bofd itaiicised terms and othsr i~;l;rr.rn:ation-de~ndent details.

1

he routme n a m mir, be precision-dependent.

~ !%!ease read the appropriate implementation dc <urnent to check the interpretation of

1. Burpose

D02BDF integrates a systsn oi' first-order ordinary differentlul eqiialíons over a

initial conditions, using a Xunge-Kutta-Merson method, and computes a gioba

stiffness check is also available.

2. Specification

SUBROUTINE D928i:F ( 8 , XEieD, N , Y , TOL,

IRELAB.

:&IJ,

STIFF,

1 YNORM, W , I % , 8 , OaTPUT, IFAIL)

I

--

I , I I Unchanged

on

exit.

(78)

DOZBDF l X 2

-

Ordinary Differential Equations

N - INTEGER.

On entry,

N

must spwcjfy the number of differential equations.

Unchanged on exit.

Y

-

real array of dimension! at least (N).

Before entry,

Y( I),Y(2),

...,

Y(N) must contain the initial values of the s

On exit,

Y(l),Y(2),

...$

computed values of the value of

T.

TOL

-

real.

On entry,

TOL

must spsi:i3j. a positive

tolerance for controlling the error in the integration. The routine D02BDF has been designed so that, for m a t problems, a

preference

in

the choice of emw test then

Page 2

IRELAB shoaild be given the value O which will result in a mixed error test.

IRELAB is unchanged on exit.

FCN

-

SUBROUT’INE, supplied by the user. FCN must evaluaie the functions Fi (i.e. the derivatives .Yí”) for given values of its

arguments

T,Y,

,...,

Y,.

Its specification

is;

S U B R O U T I N E

F C N ( T , Y . F )

real T I Y (n

1,

F (n)

where n is the actual value of

N

in the ca3 of

D02BDF.

T

-

r e d .

On

entry, T specifi

argument

T.

e of the

lts-value ;nust not

b

Y

-

r e d array of DIM

F

-

real array of DJMENSION

(n).

stiff equatícms sxch as D02EAF or DO2QBF rather than Merson’s method. If the averaged

stiffness value ar.

T

=

XEND is

much

smaller

than the values cbtained during the integration then the problem may not be stiff. Further chech should be made, possibly by varying

(79)

(80)

WZBDF

6.

Error

indicators and

N'rarnings

Errors detected by the rout IFAIL = 1 On entry, N I O or or or or STIFF. IFAIL = 2

With the given value progress can be made range from the CU Section

li).

IFAIL = 3 The dependence of th be lost if further integration range Section li). IFAIL = 4

TOL is too small for t initial integration step. retain their initial values.

Too many integratio continue to compute Section li). IFAIL = 6 The range

X

small change change in the

be treated as a warning that the computed solution is likely to be more accurate than wouid be produced by us

of TOL on a longer range error exit may occur at range is short and t specified by TOL is re large).

IFAIL = 7

An auxiliary solution used in computing the global error estimate has a component whose modulus is greater than YNORM (see Section

11).

Page 4

DO2

-

Ordifiary Differential Equations

IFAIL = 8

A serious error has occurred in an internal call to D02PAF. Check ail subroutine calls and array dimensions. Seek expert help. IFAIL = 9

The computed solution has a component whose modulus its greater than YNORM (see Section i 1)

IFAIL = 12 to 18

When STIFF

>

0.0 two solutions of the problem are computed with different tolerances (see Section 10) by internal calls of D02PAF. The error indicators IFAIL =

12,13,

....

18 conespond to the second computation and have the same meaning as IFAIL = 2,3,

...,

8 respectively. IFAIL = 16 should not occur, but if it does, the integration range should be lengthened or TOL reduced, and the computation should be restarted.

ted to sites in

8.

Timing

This depeiids on the complexity and mathematical properties of the system of differential equations defined by FCN, on the range, the tolerance

and

the value of

M.

There is an overhead of the form Ai-BXN where

A

and

1B are machine-dependent computing times. Use of this routine is likely to cost up to 50% more computing time thzac the routines D02BAF (or

D02PAF) to mìïpiite resuits of the same accuracy (though only DO2BDF gives a global error estimate). The skiffness check options given by STIFF

>

0.0 will also use more computing time; a call to D02BDF with STIFF

>

0.0 will use about twice as much time as a call to

D02BAF to comprire a result of the same accuracy.

8.

Storage

The storage required by internally declared arrays is 28 r e d elements.

10. Accuracy

The accuracy depends on TOL, on IRELAB, on the matheinatical properties of the differential system, on the length of the range of integration

and on the method. It can be controlled by

varying TOE but the approximate proportionality of the error to TOL holds only for a restricted range of values of TOL. For TOL

(81)

002

-

Ordinary Differential Equations DOLBDF

i

too large, the underlying theory may break down and the result of varying ’TOL may be unpredictable. For TOL too small, rounding error may affect the solution significantly and an error exit with IFAIL =- 2, 3, 4, 12, 13 or 14 is

possible. Internally D02BDF makes two calls to D02PAF, the first with 16.0 X TOL as error tolerance, which can lead to error exits with IFAIL

>

10.

The accuracy of the saiution computed by DO2BDF with error tolerance TOL is likely to be close to the accuracy obtainable from ~ 1 . call to

D02BAF or D02PAF with the same toierance. For small values of TOL, the global error estimate is usually accurate to within an order of magnitude. If the problem being integrated is stiff, then the theory used in designing the global error estimate breaks down and the estimaie may

be inaccurate. See the first example of Section 13

for a stiff problem where the global error is severely underestimated, When the system of equations is stiff with rapidly âecaying components, or has highly ocscillatory components, the error estimate is often an overestimate of the true error. It

is

not sufficient to make just one call to D02BDF in order to estimate the value of

TOL

which is required for a accuracy. A check should be made by

D02BDF with tire vdxies of TOL and the response of the g%cbal error estimate variation in TOL.

11. Further Comments

The routine D02BDF is not intended to be used as a general purpose integrator for differential equations. The user s

D02BAF, DQ2BBF o usually prove consid

computing time. DO2BDP should be used bo

investigate the behaviour of the error with variations in the a

perhaps as a prelim systems of different

mathematical propert course, D02BDF should also be used er a global error estimate is essential.

If all th.2t is rtrr;Uire#

i

check of stihess,

then D02BDF should d with STIFF

<

0.0, but care should b ed in interpreting the results. A more cautkais stiffness ckiack (and a global error estimate vided by a call with STIFF

>

O.O. In this he Pize of the value returned by STIFF is sure of the stiffness of the system of differential equations (see Section 5). A value greater than 0.75 should be treated as an indication that Gear’s method D02EAF or DO2QBF would

be

better than the Merson method D02BAF or DO2PAF. For smaller values of STIFF (say, STIFF

>

0.5)

Gear’s method may still be less expensive in computing time than the Merson method. If the routine fails with IFAIL = 4 or 14, then

it can be called again with a larger value of TOL if this has not already been tried. If the accuracy requested is really needed and cannot be obtained wit5 this routine, the system may be stiff or ‘oadiy scaled and possibly cannot be solved to the required accuracy.

If the routine fails with IFAIL = 2,3,

12

or 13, it is likely that DO2BDF has been called with a value for TQL which is so small that a solution cannot be obtained on the range

X

to XEND. This can happen for well-behaved systems and very small yalucs of TOL. The user should, however, consider whether there is a more fundamental difficulty:

(i) in the region of a singularity (infinite value) of the solution, the routine will normally stop with IFAIL = 2,3, 12 or 13, unless overfiow occurs first. The parameter YNORM can be used to trap the increasing solution before overflow occurs. This will lead to an exit with IFAIL = 7, 9 or 17. An exit with IFAIL

--

7 or 17 should be treated with some suspicion, as the exit with IFAiL = 9 would normally be expected to occur first. It is possible for ihese error exits to occur as a resuit of inherent instability in the differential system (see [3]). If this is suspected then a check can

be

made by changing TOL. If the same error exit occurs at appnoximately the same value of

X,

it is likely that the differential system has a singdarity. Note that numerical integration cannot be contlfiued throilgh a singularity, and analytic treatment sholald

be

considered.

(ii) the equ;tions may be stiff. This should

be evident oil a call with STIFFfO.0. If

STIFF = 0.0 on entry, then stiffness may exhibit itself by the computing time being excessively long or, occasionally,

?q aïì exit with IFAIL = 2, 3 ,

i

2 or 13.

For well-behaved systems with no difficulties such as stiffness or singularities, the Merson method should work well for low accuracy calculations (three or four figures). For high accuracy calculations,, or where FCN is costly to evaluate, Merson’s method may not be appropriate and a computationally less expensive method mag- be the Adams method D02CAF or DO2QAF. In particular an exit with IFAIL = 5 or IFAIL = 15 when the system is not stiff should be takcv as an indication that an Adams method should be used.

(82)

D02BDF ,702

- Ordinary

Differential Equations

12. Keywords

Global Error Estimate, Initial Value Problems,

Ordinary Differential Equations, Runge-Kutta--Mer it>n Method,

Stiffness Check.

_-

13.

Example

We integrate two differen

(i)

We integrate the s

y, - 0 . 0 4 ~ ~

with initial values

1

.OE-7. We print

y2 = y3 = O, across the range [O.@ 0.31 with

TOL

= 1.OE-6 and

al

error rstimate'and stiffness factor every '75 steps. Note that the is less than the intermediate stiffness factors, as is o€ten the case.

If

the e knger, the two values would be cl

iff) projectile problem

across the range [O.O, 8.û] with initial values y = O,

v

= 0.5,

0

= aj5. We compare the

solutions for

TOL

= i .OE4 and

1

.OE-5. We output the solutiûti at each integration step and include the stiffness check for comparison with

(i)

above. The value of P is obtained by using

XO

1

AAF.

13.1.

Program Text

WARNING: This single preciskm example program may require amendment for certain implementations. The results produced may not be the same. If in doubt, please seek further advice (see Essential Introductiom to the Library Manual).

DOEBDF EXAMPL

..

LOCAL ARRA REAL X O I A A F

. .

SUBROUTINE DOEBDF, OUT EXTERNAL F C N I , FCN, ObT DATA NOUT / ô / WRITE (NOUT,SS N - 3 M = 75 WRITE [NOUT,9993S) PI = XD'lAAF(X) YNORM = O . N2 = 14

DO

E O

I = 1 , 2

. .

(83)

Do2 -

Ordinary Differential Eqt.x4tiom DO2,BDF

‘-* i

TOL = 10.**(-5-1] WRITE (NOUT.93998) TOL STIFF

-

1. IR = O

x

=

o .

XEND = 0 . 3 Y ( l ) = 1 . Y ( 2 ) 5

o .

Y(3) =

o .

IFAIL = 1

CALL D02BDF(X, XEND, N, Y , TOL, IR, FCN, STIFF, YNORM. W , WRITE (NOUi,99999) IFWIL

*

N2, M . OUT, IFAIL) 20 CONTINUE WRITE (NOUT,99395] M = l DO 40 I = 1 , 2 TOL 10.**[-3-I] WRITE (NOUT,99398) TOL STIFF = 1. IR = O

x

-

o .

XEND = 8 . Y(1) =

o .

Y ( 2 ) 6 0 . 5 Y(3) = 0 . 2 ” p H IFAIL = 1

CALL û026oF(X. XFFBD, N , Y , TOL, IR, FCNI, STIFF, YNORM, W , WRITE (NOU?,99999) ZFAIL

CALL OUT(X. V . Id, STIFF)

*

NE, M . O U T , IFAIL) 40 CONTINUE

STOP

99999 FORMAT (8HO 1 , 6 8 . l ~ = , 11)

99998 FORMAT CIH0/22H CALCULATIQN WITH T W a , E g . ? )

99997 FORMAT (4(1X/), 99996 FORMAT [IX/21iiO

99995 FORMAT (1X/19HO ;ECTPLE PROBLEM)

il QQEBDF EXAMPLE PROGRAM RE$&TSEIX)

F F E

STIFF PROBLEM! END SUBROUTINE FCNáT, Y , F ) C

c

REAL

. .

F(3), Y [ Y ) F ( 1 ) = - 0 . 0 4 * L [ ? j + ?.E4*Y(2)*Y(3)

F ( 2 ) = 0.04*yl?f

-

f.E4*Y(2)*Y(3)

-

3.E7*Y(2)*F[2]

F ( 3 ) = 3.E7*Y<Zj*V(S) RETURN END SUBROUTINE FCrJ’i(7,

Y ,

F ) REAL T C

. .

SCALAR ARGUME%iS

. .

C

. .

A R R A Y ARGUMC.1TS

. .

-

INAGFLIB:I566/0: Mkl0:I 9th Jwt:ary i9831

i i

,

(84)

DûLBDF Dû.2

-

Ordinary Differential Equations REAL F ( 3 ) , Y ( 3 ) C

. .

C

. .

FUNCTION REFERENCES

. .

C

. .

REAL COS, S I N F(1) = S I N ( Y ( 3 ) ) / C O S ( Y ( 3 ) ) F ( 2 ) -O.O32*F(1)/ì'(2)

-

O.O2*Y(2)/COS(Y(3)) RETURN END SUBROUTINE OUT(X, Y , W , S T l F F ) REAL S T I F F , X REAL W(3,14). Y[3) F ( 3 ) 5 -O.O32/(Y 2 )

1

C

, .

SCALAR ARGUMENTS

. .

C

. .

ARRAY ARGUME C

. .

C

. .

LOCAL SCALA INTEGER J, NOU C

. .

DATA WOUT /6/

WRITE (NOUT,99999) X , ( Y ( J ) , J = 1 . 3 ) , (W(J,'I),J=1,3), RETURN

*

(W(J,S), J=1,3), [WlJ, 3 ) , J = l , 3 ) , S T I F F

99999 FORMAT (15HOX AND SOLUTION, F 1 3 . 5 , 3E13.5116H CURRENT ERROR E ,

*

IIHSTIMATES , 3E12.4127H MAXIMUM ERROR ESTIMATES

*

3E12.4/41H NUMBER OF SIGN CHANGES FOR EACH ESTIMATE,

*

3F5.0/17H STIFFWES? FACTOR, 3E11.4)

END

13.2. Program

Data

None.

13.3. Program Results

DO2BDF EXAMPLE PRO

SIMPLE S T I F F PROBLE CALCULATION WITH T

X AND SOLUTION .I3795 0.99463E+00 0.35541E-04 0.53355E-02 CURRENT ERROR E S T I 0.5721E-06 0.1407E-05 -0.1980E-05 MAXIMUM ERROR ESTIUAT 0.6003E-06 0.1886E-05 -0.2995E-05

NUgBER ^F SIV3

Y H A

i. 7 î . 35.

STIFFNESS FACTOR O

X AND SOLUTION 0.27370 O.Y8962E+00 0.34644E-G4 C.10349E-01 CURRENT ERROR E S T I 0.2037E-05 -0.2890E-36 -5.1748E-05 MAXIMUM ERROR E S T I 0.2037E-05 0.1886E-05 -0.3282E-05 NUMBER OF

SIGN

CHANGES FOR EACH ESTIMATE 1 . 148. 35.

STIFFNESS FACTOR 0.10r30E+01 IFAIL=O

X AND SOLUTION 0.30000 0.98867E+00 0.34477E-O4 0.11292E-01 CURRENT ERROR ESTIgATES 0.2395E-05 -0.2408E-07 -0.237OE-O5 MAXIMUM ERROR ESTIMATES 0.2406E-05 0.1886E-05 -0.3693E-05 NUMBER OF

SIGN

CHANGES FOR EACH ESTIMATE I. 160. 35.

(85)

Do2

-

Ordinary Differentia! Cqudons

~

S T I F F N E S S FACTOR 0.8967L+00 CALCULATION WITH TOL= G.1E-06

X AND SOLUTION 0.72268 0.99521E+00 0.35647E-U4 0.47550E-02 CURRENT ERROR ESTIMATES 0.3712E-08 -0.8869E-08 0.5157E-08 MAXIMUM ERROR ESTIMATES 0.5679E-08 -0.1488E-06 -0.1475E-O6 NUMBER OF S I G N CHANGES FOR EACH ESTIMATE 19. 7 2 . 72.

S T I F F N E S S FACTOR C.?000E+01

X AND SOLUTION 0.24BOE 0.99062E+00 0.34822E-04 0.93475E-02 CURRENT ERROR ESTIMATES 0.1196E-07 0.1784E-07 -0.2980E-07 MAXIMUM ERROR ESTSRATES 0.1536E-07 0.1535E-06 -0.1558E-06 NUMBER OF S I G N C H A N E S FOR EACH ESTIMATE 19. 147. 147. S T I F F N E S S FACTOR û.'iOOL'1E+OI

I F A I L m O

X AND SOLUTION 0.30000 0.98867E+00 0.34477E-04 0.11292E-01 CURRENT ERROR ESTIMATES 0.1880E-07 0.5409E-08 -0.2421E-07 MAXIMUM ERROR ESTIMATES 0.2456E-07 0.1535E-06 -0.1358E-O6 NUMBER OF S I G N CHANGES FOR EACH ESTIMATE 1 9 . 179. 177.

S T I F F N E S S FACTOR 0.989ilE+00

PROJECTILE PROBLEM

CALCULATION WITH . 'IE-03

CURRENT ERROR EST O 1146E-03 -0.6312E-05 -0.1450E-04 MAXIMUM ERROR ESTIMATES 0.1146E-03 -0.6312E-05 -0.3450E-04 NUMBER OF SIGN CHANGES FOR EACH ESTIMATE - 1 . -1. -1.

X AND SOLUTION 07212 0.96785E+00 0.38894E+00 -0.36672E+00 CURRENT ERROR ESTIMATES 0.4887E-O4 0.4311E-06 -0.L230E-04 MAXIMUM ERROR ESTIMATES O . 1146E-03 -0.6312E-05 -0.3230E-04 NUMBER OF SIGN CHARGES FOR EACH ESTIMATE O . 1 . - 4 I .

S T I F F N E S S FACTOR

CURRENT ERROR E S -0.9947E-04 0.8658E-05 -0.3092E-04

UAliii#üY ERROR

ES

U.1146E-03 0.8658E-05 -0.3230E-04

NUMBER OF S I G N CHAMGES FOR EACH ESTIMATE 1 . 1 . O .

X AND SOLUTION 8.011000 -0.12458E+01 0.51298E+00 -0.85368E+00 CURRENT ERROR ESTIMATES -0.1535E-03 0.1035E-O4 -0.2605E-04 MAXIMUM ERROR E S -0.1535E-03 0.1035E-04 -0.3230E-04 NUMBER OF S I G N C EACH ESTIMATE 1 . I . O .

STIFFNESS FACTOR

X AND SOLUTION 8.00000 -0.12458E+01 0.51298E+00 -0.85368E+00 CURRENT ERROR ESTIMATES -0.1535E-03 0.1035E-04 -0.2605E-04 MAXIMUM ERROR ESTIMATES -0.1535E-03 0.1035E-04 -0.3230E-04 NUMBER OF S I G N C EACH ESTIMATE 1 . 1 . O .

CALCULATION WITH TOL= 0 . 1 E - 0 4

X AND SOLUTION 3.54294 0.8824GE+00 0.42241L+00 0.39212E+00 CURRENT ERROR ESTIMATES 0.4705E-05 -0.2472E-06 -0.5015E-06 MAXIMUM ERROR ESTIMATES 0.4705E-05 -0.2472E-06 -0.5013E-06 NUMBER OF S I G N CrlANGES

F O R

E k f H ESTIMATE - 1 . - 1 . - 1 .

X AND SOLUTION 4337 0.12248E+01 0.385255+00 0.15374E+00

S T I F F N E S S FACTOR 00E+Ilil X AND SOLUTION 07048 - o . z 9 9 2 2 ~ + 0 0 0 . 4 6 3 ~ 5 ~ + 3 0 - 0 . 7 2 8 4 3 ~ + 0 0 S T I F F N E S S FACTOR 0.0000E+00 I F A I L = O L.- [NAGFLìB:1566/0:MklO:l~t~~ Januoy 19831 D02BDF Page 9 I

(86)

DûLBDF DO2 -- Ordinary Differential Equations

S T I F F N E S S FACTOR 0.0000E+00

X AND SOLUTION 2.79376 0.12324E+01 0.38423E+00 0.14283E+00 CURRENT ERROR E S T I M A X S 0.7756E-05 -0.4082E-06 -0.1192E-05 MAXIMUM ERROR ESTIMATES 0.7756E-05 -0.4082E-06 -0.1192E-05 NUMBER OF SIGN CHANGES

F O R

EACH ESTIMATE -1. - 1 . - 1 .

X AND SOLUTION 4 . Q 0 3 3 5 0.12419E+01 0.37434E+00 -0.12968E+00 CURRENT ERROR ESTIMATES 0.7586E-05 -0.3134E-06 -0.1907E-05 MAXIMUM ERROR ESTIMATES 0.7756E-05 -0.4082E-06 -0.1907E-05 NUMBER OF SIGN CHANG EACH ESTIMATE O . O . - 1 .

S T I F F N E S S FACTOR 0 . 0

X AND SOLUTION 4 3 !3.92937E+00 0.39129E+QO -0.38722E+00 CURRENT ERROR ESTIMA 0.2269E-05 0.1091E-06 -0.2408E-05 MAXIMUM ERROR ESTIMA 0.7756E-05 -0.4082E-06 -0.2408E-05 NUMBER OF S I G N CHANGES FOR EACH ESTIMATE O . 1 . -’i.

S T I F F N E S S FACTOR O .

X AND SOLUTION 5864 0.33822E+00 0.42738E+00 -0.53711E+00 CURRENT ERROR ESTIMA -0.3846E-05 0.4791E-06 -0.2392E-05 MAXIMUM ERROR ESTIMATES 0.7756E-05 0.4791E-06 -0.2408E-05 NUMBER OF SIGN CHANGES FOR EACH ESTIMATE 1 . 1 . O .

S T I F F N E S S FACTOR 0.00QGE+O0

X AND SOLUTION -0.54587E+00 0.47731E+@O -0.76687E+00 CURRENT ERROR E S T I M -0.9336E-05 0.7032E-06 -0.2068E-05 MAXiMUM ERROR E S T I M -0.9336E-05 0.7032E-06 -0.2408E-O5 NUMBER OF SIGN CHAN EACH ESTIMATE 1 . 1 . O .

O . O B U O E + O G

X AND SOLUTION O0 -0.12459E+01 0.51299E+00 -0.85371E+00 CURRENT ERROR E S T I M -0.1208E-04 0.7770E-06 -0.1810E-05 MAXIMUM ERROR E S T I M -0.1208E-04 0.7770E-06 -0.2403E-05 NUMBER OF SIGN CHANGES FOR EACH ESTIMATE 1 . 1 . O .

S T I F F N E S S FACTOR O .

X AND SOLUTION 00130 -0.12459E+01 0.51299E+OO ;0.85371E+00

CURRENT ERROR E S T I M -0.1208E-04 0.7770E-06 -0.18’lOE:05 MAXIMUM ERROR E S T I M -0.1208E-04 0.7770E-06 -0.2408E-05 NUMBER

OF

SIGN CHAN EACH ESTIMATE 1 . 1 . O .

O.OOOOE+OO

S T I F F N E S S FACTOR 0.0000€+00

IFAIL=O