Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DormandPrinceRK56 Class Reference

G4DormandPrinceRK56 implements the 6(5) embedded Runge-Kutta non-FSAL method. More...

#include <G4DormandPrinceRK56.hh>

Inheritance diagram for G4DormandPrinceRK56:

Public Member Functions

 G4DormandPrinceRK56 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 ~G4DormandPrinceRK56 () override
 G4DormandPrinceRK56 (const G4DormandPrinceRK56 &)=delete
G4DormandPrinceRK56 & operator= (const G4DormandPrinceRK56 &)=delete
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
G4double DistChord () const override
G4int IntegratorOrder () const override
G4StepperType StepperType () const override
void SetupInterpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step)
void SetupInterpolate (const G4double yInput[], const G4double dydx[], const G4double Step)
void SetupInterpolation ()
void Interpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
void Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
void Interpolate (G4double tau, G4double yOut[])
void SetupInterpolate_high (const G4double yInput[], const G4double dydx[], const G4double Step)
void Interpolate_high (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
virtual ~G4MagIntegratorStepper ()=default
 G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
G4MagIntegratorStepper & operator= (const G4MagIntegratorStepper &)=delete
void NormaliseTangentVector (G4double vec[6])
void NormalisePolarizationVector (G4double vec[12])
void RightHandSide (const G4double y[], G4double dydx[]) const
void RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
G4int GetNumberOfVariables () const
G4int GetNumberOfStateVariables () const
G4int IntegrationOrder ()
G4EquationOfMotion * GetEquationOfMotion ()
const G4EquationOfMotion * GetEquationOfMotion () const
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long GetfNoRHSCalls ()
void ResetfNORHSCalls ()
G4bool IsFSAL () const
G4bool isQSS () const
void SetIsQSS (G4bool val)

Additional Inherited Members

Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (G4int order)
void SetFSAL (G4bool flag=true)

Detailed Description

G4DormandPrinceRK56 implements the 6(5) embedded Runge-Kutta non-FSAL method.

Definition at line 45 of file G4DormandPrinceRK56.hh.

Constructor & Destructor Documentation

◆ G4DormandPrinceRK56() [1/2]

G4DormandPrinceRK56::G4DormandPrinceRK56 ( G4EquationOfMotion * EqRhs,
G4int numberOfVariables = 6,
G4bool primary = true )

Constructor for G4DormandPrinceRK56.

Parameters
[in]EqRhsPointer to the provided equation of motion.
[in]numberOfVariablesThe number of integration variables.
[in]primaryFlag for initialisation of the auxiliary stepper.

Definition at line 37 of file G4DormandPrinceRK56.cc.

40 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
41{
42 const G4int numberOfVariables = noIntegrationVariables;
43
44 // New Chunk of memory being created for use by the stepper
45
46 // aki - for storing intermediate RHS
47 //
48 ak2 = new G4double[numberOfVariables];
49 ak3 = new G4double[numberOfVariables];
50 ak4 = new G4double[numberOfVariables];
51 ak5 = new G4double[numberOfVariables];
52 ak6 = new G4double[numberOfVariables];
53 ak7 = new G4double[numberOfVariables];
54 ak8 = new G4double[numberOfVariables];
55 ak9 = new G4double[numberOfVariables];
56
57 // Memory for Additional stages
58 //
59 ak10 = new G4double[numberOfVariables];
60 ak11 = new G4double[numberOfVariables];
61 ak12 = new G4double[numberOfVariables];
62 ak10_low = new G4double[numberOfVariables];
63
64 const G4int numStateVars = std::max(noIntegrationVariables, 8);
65 yTemp = new G4double[numStateVars];
66 yIn = new G4double[numStateVars] ;
67
68 fLastInitialVector = new G4double[numStateVars] ;
69 fLastFinalVector = new G4double[numStateVars] ;
70
71 fLastDyDx = new G4double[numStateVars];
72
73 fMidVector = new G4double[numStateVars];
74 fMidError = new G4double[numStateVars];
75
76 if( primary )
77 {
78 fAuxStepper = new G4DormandPrinceRK56(EqRhs, numberOfVariables, !primary);
79 }
80}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4DormandPrinceRK56::G4DormandPrinceRK56
G4DormandPrinceRK56(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition G4DormandPrinceRK56.cc:37
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
Definition G4MagIntegratorStepper.cc:37

Referenced by G4DormandPrinceRK56(), G4DormandPrinceRK56(), and operator=().

◆ ~G4DormandPrinceRK56()

G4DormandPrinceRK56::~G4DormandPrinceRK56 ( )
override

Destructor.

Definition at line 84 of file G4DormandPrinceRK56.cc.

85{
86 // clear all previously allocated memory for stepper and DistChord
87
88 delete [] ak2;
89 delete [] ak3;
90 delete [] ak4;
91 delete [] ak5;
92 delete [] ak6;
93 delete [] ak7;
94 delete [] ak8;
95 delete [] ak9;
96
97 delete [] ak10;
98 delete [] ak10_low;
99 delete [] ak11;
100 delete [] ak12;
101
102 delete [] yTemp;
103 delete [] yIn;
104
105 delete [] fLastInitialVector;
106 delete [] fLastFinalVector;
107 delete [] fLastDyDx;
108 delete [] fMidVector;
109 delete [] fMidError;
110
111 delete fAuxStepper;
112}

◆ G4DormandPrinceRK56() [2/2]

G4DormandPrinceRK56::G4DormandPrinceRK56 ( const G4DormandPrinceRK56 & )
delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

◆ DistChord()

G4double G4DormandPrinceRK56::DistChord ( ) const
overridevirtual

Returns the distance from chord line.

Implements G4MagIntegratorStepper.

Definition at line 357 of file G4DormandPrinceRK56.cc.

358{
359 G4double distLine, distChord;
360 G4ThreeVector initialPoint, finalPoint, midPoint;
361
362 // Store last initial and final points
363 // (they will be overwritten in self-Stepper call!)
364 //
365 initialPoint = G4ThreeVector( fLastInitialVector[0],
366 fLastInitialVector[1], fLastInitialVector[2]);
367 finalPoint = G4ThreeVector( fLastFinalVector[0],
368 fLastFinalVector[1], fLastFinalVector[2]);
369
370 // Do half a step using StepNoErr
371
372 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
373 fMidVector, fMidError );
374
375 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
376
377 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
378 // distance of Chord
379 //
380 if (initialPoint != finalPoint)
381 {
382 distLine = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
383 distChord = distLine;
384 }
385 else
386 {
387 distChord = (midPoint-initialPoint).mag();
388 }
389 return distChord;
390}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition G4LineSection.hh:96

◆ IntegratorOrder()

G4int G4DormandPrinceRK56::IntegratorOrder ( ) const
inlineoverridevirtual

Returns the order, 5, of integration.

Implements G4MagIntegratorStepper.

Definition at line 95 of file G4DormandPrinceRK56.hh.

95{ return 5; }

◆ Interpolate() [1/2]

void G4DormandPrinceRK56::Interpolate ( const G4double yInput[],
const G4double dydx[],
const G4double Step,
G4double yOut[],
G4double tau )
inline

Wrappers for Interpolate_low() above.

Definition at line 136 of file G4DormandPrinceRK56.hh.

141 {
142 Interpolate_low( yInput, dydx, Step, yOut, tau);
143 }
G4DormandPrinceRK56::Interpolate_low
void Interpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition G4DormandPrinceRK56.cc:432

Referenced by Interpolate().

◆ Interpolate() [2/2]

void G4DormandPrinceRK56::Interpolate ( G4double tau,
G4double yOut[] )
inline

Definition at line 144 of file G4DormandPrinceRK56.hh.

145 {
146 Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau );
147 }
G4DormandPrinceRK56::Interpolate
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition G4DormandPrinceRK56.hh:136

◆ Interpolate_high()

void G4DormandPrinceRK56::Interpolate_high ( const G4double yInput[],
const G4double dydx[],
const G4double Step,
G4double yOut[],
G4double tau )

Calculates the output at the tau fraction of Step, using the polynomial coefficients and the respective stages.

Definition at line 569 of file G4DormandPrinceRK56.cc.

574{
575 // Define the coefficients for the polynomials
576 //
577 G4double bi[13][6], b[13];
578 G4int numberOfVariables = GetNumberOfVariables();
579
580
581 // COEFFICIENTS OF bi[ 1]
582 bi[1][0] = 1.0 ,
583 bi[1][1] = -18487.0/2880.0 ,
584 bi[1][2] = 139189.0/7200.0 ,
585 bi[1][3] = -53923.0/1800.0 ,
586 bi[1][4] = 13811.0/600,
587 bi[1][5] = -2071.0/300,
588 // --------------------------------------------------------
589 //
590 // COEFFICIENTS OF bi[2]
591 bi[2][0] = 0.0 ,
592 bi[2][1] = 0.0 ,
593 bi[2][2] = 0.0 ,
594 bi[2][3] = 0.0 ,
595 bi[2][4] = 0.0 ,
596 bi[2][5] = 0.0 ,
597 // --------------------------------------------------------
598 //
599 // COEFFICIENTS OF bi[3]
600 bi[3][0] = 0.0 ,
601 bi[3][1] = 0.0 ,
602 bi[3][2] = 0.0 ,
603 bi[3][3] = 0.0 ,
604 bi[3][4] = 0.0 ,
605 bi[3][5] = 0.0 ,
606 // --------------------------------------------------------
607 //
608 // COEFFICIENTS OF bi[4]
609 bi[4][0] = 0.0 ,
610 bi[4][1] = -30208.0/14157.0 ,
611 bi[4][2] = 1147904.0/70785.0 ,
612 bi[4][3] = -241664.0/5445.0 ,
613 bi[4][4] = 241664.0/4719.0 ,
614 bi[4][5] = -483328.0/23595.0 ,
615 // --------------------------------------------------------
616 //
617 // COEFFICIENTS OF bi[5]
618 bi[5][0] = 0.0 ,
619 bi[5][1] = -177147.0/32912.0 ,
620 bi[5][2] = 3365793.0/82280.0 ,
621 bi[5][3] = -2302911.0/20570.0 ,
622 bi[5][4] = 531441.0/4114.0 ,
623 bi[5][5] = -531441.0/10285.0 ,
624 // --------------------------------------------------------
625 //
626 // COEFFICIENTS OF bi[6]
627 bi[6][0] = 0.0 ,
628 bi[6][1] = 536.0/141.0 ,
629 bi[6][2] = -20368.0/705.0 ,
630 bi[6][3] = 55744.0/705.0 ,
631 bi[6][4] = -4288.0/47.0 ,
632 bi[6][5] = 8576.0/235,
633 // --------------------------------------------------------
634 //
635 // COEFFICIENTS OF bi[7]
636 bi[7][0] = 0.0 ,
637 bi[7][1] = -1977326743.0/723932352.0 ,
638 bi[7][2] = 37569208117.0/1809830880.0 ,
639 bi[7][3] = -1977326743.0/34804440.0 ,
640 bi[7][4] = 1977326743.0/30163848.0 ,
641 bi[7][5] = -1977326743.0/75409620.0 ,
642 // --------------------------------------------------------
643 //
644 // COEFFICIENTS OF bi[8]
645 bi[8][0] = 0.0 ,
646 bi[8][1] = 259.0/144.0 ,
647 bi[8][2] = -4921.0/360.0 ,
648 bi[8][3] = 3367.0/90.0 ,
649 bi[8][4] = -259.0/6.0 ,
650 bi[8][5] = 259.0/15.0 ,
651 // --------------------------------------------------------
652 //
653 // COEFFICIENTS OF bi[9]
654 bi[9][0] = 0.0 ,
655 bi[9][1] = 62.0/105.0 ,
656 bi[9][2] = -2381.0/525.0 ,
657 bi[9][3] = 949.0/75.0 ,
658 bi[9][4] = -2636.0/175.0 ,
659 bi[9][5] = 1112.0/175.0 ,
660 // --------------------------------------------------------
661 //
662 // COEFFICIENTS OF bi[10]
663 bi[10][0] = 0.0 ,
664 bi[10][1] = 43.0/3.0 ,
665 bi[10][2] = -1534.0/15.0 ,
666 bi[10][3] = 3767.0/15.0 ,
667 bi[10][4] = -1264.0/5.0 ,
668 bi[10][5] = 448.0/5.0 ,
669 // --------------------------------------------------------
670 //
671 // COEFFICIENTS OF bi[11]
672 bi[11][0] = 0.0 ,
673 bi[11][1] = 63.0/5.0 ,
674 bi[11][2] = -1494.0/25.0 ,
675 bi[11][3] = 2907.0/25.0 ,
676 bi[11][4] = -2592.0/25.0 ,
677 bi[11][5] = 864.0/25.0 ,
678 // --------------------------------------------------------
679 //
680 // COEFFICIENTS OF bi[12]
681 bi[12][0] = 0.0 ,
682 bi[12][1] = -576.0/35.0 ,
683 bi[12][2] = 19584.0/175.0 ,
684 bi[12][3] = -6336.0/25.0 ,
685 bi[12][4] = 41472.0/175.0 ,
686 bi[12][5] = -13824.0/175.0 ;
687 // --------------------------------------------------------
688
689 for(G4int i = 0; i< numberOfVariables; ++i)
690 {
691 yIn[i] = yInput[i];
692 }
693
694 G4double tau0 = tau;
695
696 // Calculating the polynomials (coefficents for the respective stages) :
697 //
698 for(auto i=1; i<=12; ++i) // i is NOT the coordinate no., it's stage no.
699 {
700 b[i] = 0;
701 tau = 1.0;
702 for(auto j=0; j<=5; ++j)
703 {
704 b[i] += bi[i][j]*tau;
705 tau*=tau0;
706 }
707 }
708
709 // Calculating the interpolation at the fraction tau of the step using
710 // the polynomial coefficients and the respective stages
711 //
712 for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
713 {
714 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
715 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
716 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
717 b[10]*ak10[i] + b[11]*ak11[i] + b[12]*ak12[i]);
718 }
719}
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

◆ Interpolate_low()

void G4DormandPrinceRK56::Interpolate_low ( const G4double yInput[],
const G4double dydx[],
const G4double Step,
G4double yOut[],
G4double tau )

Calculates the output at the tau fraction of Step.

Definition at line 432 of file G4DormandPrinceRK56.cc.

437{
438 G4double bf1, bf4, bf5, bf6, bf7, bf8, bf9, bf10;
439
440 G4double tau0 = tau;
441 const G4int numberOfVariables= this->GetNumberOfVariables();
442
443 for(G4int i=0; i<numberOfVariables; ++i)
444 {
445 yIn[i]=yInput[i];
446 }
447
448 G4double tau_2 = tau0*tau0 ,
449 tau_3 = tau0*tau_2,
450 tau_4 = tau_2*tau_2;
451
452 // bf2 = bf3 = 0.0
453 bf1 = (66480.0*tau_4-206243.0*tau_3+237786.0*tau_2-124793.0*tau+28800.0)
454 / 28800.0 ;
455 bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*tau_2 + 119706.0*tau -40973.0)
456 / 70785.0 ;
457 bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*tau_2 + 34786.0*tau - 9293.0)
458 / 1645600.0 ;
459 bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 23786.0*tau - 6533.0)
460 / 705.0 ;
461 bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797.0*tau_2 + 17574.0*tau - 1927.0)
462 / 7239323520.0 ;
463 bf8 = 37.0*tau*(336.0*tau_3 - 661.0*tau_2 + 342.0*tau -31.0)
464 / 1440.0 ;
465 bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau +3.0)
466 / 4.0 ;
467 bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0*tau - 1.0) ;
468
469 for( G4int i=0; i<numberOfVariables; ++i)
470 {
471 yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i] + bf4*ak4[i] + bf5*ak5[i] +
472 bf6*ak6[i] + bf7*ak7[i] + bf8*ak8[i] +
473 bf9*ak9[i] + bf10*ak10_low[i] ) ;
474 }
475}

Referenced by Interpolate().

◆ operator=()

G4DormandPrinceRK56 & G4DormandPrinceRK56::operator= ( const G4DormandPrinceRK56 & )
delete

◆ SetupInterpolate()

void G4DormandPrinceRK56::SetupInterpolate ( const G4double yInput[],
const G4double dydx[],
const G4double Step )
inline

Wrappers for SetupInterpolate_low() above.

Definition at line 113 of file G4DormandPrinceRK56.hh.

116 {
117 SetupInterpolate_low( yInput, dydx, Step);
118 }
G4DormandPrinceRK56::SetupInterpolate_low
void SetupInterpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition G4DormandPrinceRK56.cc:401

Referenced by SetupInterpolation().

◆ SetupInterpolate_high()

void G4DormandPrinceRK56::SetupInterpolate_high ( const G4double yInput[],
const G4double dydx[],
const G4double Step )

Prepares the interpolant and calculates the extra stages. Sixth order interpolant with 3 additional stages per step.

Definition at line 489 of file G4DormandPrinceRK56.cc.

492{
493 // Coefficients for the additional stages
494 //
495 G4double b101 = 33797.0/460800.0 ,
496 b102 = 0.0 ,
497 b103 = 0.0 ,
498 b104 = 27757.0/70785.0 ,
499 b105 = 7923501.0/26329600.0 ,
500 b106 = -927.0/3760.0 ,
501 b107 = -3314760575.0/23165835264.0 ,
502 b108 = 2479.0/23040.0 ,
503 b109 = 1.0/64.0 ,
504
505 b111 = 5843.0/76800.0 ,
506 b112 = 0.0 ,
507 b113 = 0.0 ,
508 b114 = 464.0/2673.0 ,
509 b115 = 353997.0/1196800.0 ,
510 b116 = -15068.0/57105.0 ,
511 b117 = -282475249.0/3644974080.0 ,
512 b118 = 8678831.0/156245760.0 ,
513 b119 = 116113.0/11718432.0 ,
514 b1110 = -25.0/243.0 ,
515
516 b121 = 15088049.0/199065600.0 ,
517 b122 = 0.0 ,
518 b123 = 0.0 ,
519 b124 = 2.0/5.0 ,
520 b125 = 92222037.0/268083200.0 ,
521 b126 = -433420501.0/1528586640.0 ,
522 b127 = -11549242677007.0/83630285291520.0 ,
523 b128 = 2725085557.0/26167173120.0 ,
524 b129 = 235429367.0/16354483200.0 ,
525 b1210 = -90924917.0/1040739840.0 ,
526 b1211 = -271149.0/21414400.0 ;
527
528 const G4int numberOfVariables = GetNumberOfVariables();
529
530 // Saving yInput because yInput and yOut can be aliases for same array
531 //
532 for(G4int i=0; i<numberOfVariables; ++i)
533 {
534 yIn[i]=yInput[i];
535 }
536 yTemp[7] = yIn[7];
537
538 // Evaluate the extra stages
539 //
540 for(G4int i=0; i<numberOfVariables; ++i)
541 {
542 yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
543 b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
544 b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
545 }
546 RightHandSide(yTemp, ak10); // 10th Stage
547
548 for(G4int i=0; i<numberOfVariables; ++i)
549 {
550 yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
551 b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
552 b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
553 b1110*ak10[i]);
554 }
555 RightHandSide(yTemp, ak11); // 11th Stage
556
557 for(G4int i=0; i<numberOfVariables; ++i)
558 {
559 yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
560 b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
561 b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
562 b1210*ak10[i] + b1211*ak11[i]);
563 }
564 RightHandSide(yTemp, ak12); // 12th Stage
565}
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const

◆ SetupInterpolate_low()

void G4DormandPrinceRK56::SetupInterpolate_low ( const G4double yInput[],
const G4double dydx[],
const G4double Step )

Prepares the interpolant and calculates the extra stages. Fifth order interpolant with one extra function evaluation per step.

Definition at line 401 of file G4DormandPrinceRK56.cc.

404{
405 const G4int numberOfVariables= this->GetNumberOfVariables();
406
407 G4double b_101 = 33797.0/460800.0 ,
408 b_102 = 0. ,
409 b_103 = 0. ,
410 b_104 = 27757.0/70785.0 ,
411 b_105 = 7923501.0/26329600.0 ,
412 b_106 = -927.0/3760.0 ,
413 b_107 = -3314760575.0/23165835264.0 ,
414 b_108 = 2479.0/23040.0 ,
415 b_109 = 1.0/64.0 ;
416
417 for(G4int i=0; i<numberOfVariables; ++i)
418 {
419 yIn[i]=yInput[i];
420 }
421
422
423 for(G4int i=0; i<numberOfVariables; ++i)
424 {
425 yTemp[i] = yIn[i] + Step*(b_101*dydx[i] + b_102*ak2[i] + b_103*ak3[i] +
426 b_104*ak4[i] + b_105*ak5[i] + b_106*ak6[i] +
427 b_107*ak7[i] + b_108*ak8[i] + b_109*ak9[i]);
428 }
429 RightHandSide(yTemp, ak10_low); // 10th Stage
430}

Referenced by SetupInterpolate().

◆ SetupInterpolation()

void G4DormandPrinceRK56::SetupInterpolation ( )
inline

Definition at line 119 of file G4DormandPrinceRK56.hh.

120 {
121 SetupInterpolate( fLastInitialVector, fLastDyDx, fLastStepLength);
122 }
G4DormandPrinceRK56::SetupInterpolate
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition G4DormandPrinceRK56.hh:113

◆ Stepper()

void G4DormandPrinceRK56::Stepper ( const G4double y[],
const G4double dydx[],
G4double h,
G4double yout[],
G4double yerr[] )
overridevirtual

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'h'. Integrates ODE starting values y[0 to 6]. Outputs yout[] and its estimated error yerr[].

Parameters
[in]yStarting values array of integration variables.
[in]dydxDerivatives array.
[in]hThe given step size.
[out]youtIntegration output.
[out]yerrThe estimated error.

Implements G4MagIntegratorStepper.

Definition at line 119 of file G4DormandPrinceRK56.cc.

126{
127 G4int i;
128
129 // The various constants defined on the basis of butcher tableu
130 // Old Coefficients from
131 // P.J.Prince and J.R.Dormand, "High order embedded Runge-Kutta formulae"
132 // Journal of Computational and Applied Math., vol.7, no.1, pp.67-75, 1980.
133 //
134 const G4double b21 = 1.0/10.0 ,
135 b31 = -2.0/81.0 ,
136 b32 = 20.0/81.0 ,
137
138 b41 = 615.0/1372.0 ,
139 b42 = -270.0/343.0 ,
140 b43 = 1053.0/1372.0 ,
141
142 b51 = 3243.0/5500.0 ,
143 b52 = -54.0/55.0 ,
144 b53 = 50949.0/71500.0 ,
145 b54 = 4998.0/17875.0 ,
146
147 b61 = -26492.0/37125.0 ,
148 b62 = 72.0/55.0 ,
149 b63 = 2808.0/23375.0 ,
150 b64 = -24206.0/37125.0 ,
151 b65 = 338.0/459.0 ,
152
153 b71 = 5561.0/2376.0 ,
154 b72 = -35.0/11.0 ,
155 b73 = -24117.0/31603.0 ,
156 b74 = 899983.0/200772.0 ,
157 b75 = -5225.0/1836.0 ,
158 b76 = 3925.0/4056.0 ,
159
160 b81 = 465467.0/266112.0 ,
161 b82 = -2945.0/1232.0 ,
162 b83 = -5610201.0/14158144.0 ,
163 b84 = 10513573.0/3212352.0 ,
164 b85 = -424325.0/205632.0 ,
165 b86 = 376225.0/454272.0 ,
166 b87 = 0.0 ,
167
168 c1 = 61.0/864.0 ,
169 c2 = 0.0 ,
170 c3 = 98415.0/321776.0 ,
171 c4 = 16807.0/146016.0 ,
172 c5 = 1375.0/7344.0 ,
173 c6 = 1375.0/5408.0 ,
174 c7 = -37.0/1120.0 ,
175 c8 = 1.0/10.0 ,
176
177 b91 = 61.0/864.0 ,
178 b92 = 0.0 ,
179 b93 = 98415.0/321776.0 ,
180 b94 = 16807.0/146016.0 ,
181 b95 = 1375.0/7344.0 ,
182 b96 = 1375.0/5408.0 ,
183 b97 = -37.0/1120.0 ,
184 b98 = 1.0/10.0 ,
185
186 dc1 = c1 - 821.0/10800.0 ,
187 dc2 = c2 - 0.0 ,
188 dc3 = c3 - 19683.0/71825,
189 dc4 = c4 - 175273.0/912600.0 ,
190 dc5 = c5 - 395.0/3672.0 ,
191 dc6 = c6 - 785.0/2704.0 ,
192 dc7 = c7 - 3.0/50.0 ,
193 dc8 = c8 - 0.0 ,
194 dc9 = 0.0;
195
196
197// New Coefficients obtained from
198// Table 3 RK6(5)9FM with corrected coefficients
199//
200// T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
201// "Continuous approximation with embedded Runge-Kutta methods"
202// Applied Numerical Mathematics, vol. 22, no. 1, pp. 51-62, 1996.
203//
204// b21 = 1.0/9.0 ,
205//
206// b31 = 1.0/24.0 ,
207// b32 = 1.0/8.0 ,
208//
209// b41 = 1.0/16.0 ,
210// b42 = 0.0 ,
211// b43 = 3.0/16.0 ,
212//
213// b51 = 280.0/729.0 ,
214// b52 = 0.0 ,
215// b53 = -325.0/243.0 ,
216// b54 = 1100.0/729.0 ,
217//
218// b61 = 6127.0/14680.0 ,
219// b62 = 0.0 ,
220// b63 = -1077.0/734.0 ,
221// b64 = 6494.0/4037.0 ,
222// b65 = -9477.0/161480.0 ,
223//
224// b71 = -13426273320.0/14809773769.0 ,
225// b72 = 0.0 ,
226// b73 = 4192558704.0/2115681967.0 ,
227// b74 = 14334750144.0/14809773769.0 ,
228// b75 = 117092732328.0/14809773769.0 ,
229// b76 = -361966176.0/40353607.0 ,
230//
231// b81 = -2340689.0/1901060.0 ,
232// b82 = 0.0 ,
233// b83 = 31647.0/13579.0 ,
234// b84 = 253549596.0/149518369.0 ,
235// b85 = 10559024082.0/977620105.0 ,
236// b86 = -152952.0/12173.0 ,
237// b87 = -5764801.0/186010396.0 ,
238//
239// b91 = 203.0/2880.0 ,
240// b92 = 0.0 ,
241// b93 = 0.0 ,
242// b94 = 30208.0/70785.0 ,
243// b95 = 177147.0/164560.0 ,
244// b96 = -536.0/705.0 ,
245// b97 = 1977326743.0/3619661760.0 ,
246// b98 = -259.0/720.0 ,
247//
248//
249// dc1 = 36567.0/458800.0 - b91,
250// dc2 = 0.0 - b92,
251// dc3 = 0.0 - b93,
252// dc4 = 9925984.0/27063465.0 - b94,
253// dc5 = 85382667.0/117968950.0 - b95,
254// dc6 = - 310378.0/808635.0 - b96 ,
255// dc7 = 262119736669.0/345979336560.0 - b97,
256// dc8 = - 1.0/2.0 - b98 ,
257// dc9 = -101.0/2294.0 ;
258
259 // end of declaration
260
261 const G4int numberOfVariables = GetNumberOfVariables();
262
263 // The number of variables to be integrated over
264 //
265 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
266
267 // Saving yInput because yInput and yOut can be aliases for same array
268 //
269 for(i=0; i<numberOfVariables; ++i)
270 {
271 yIn[i]=yInput[i];
272 }
273 // RightHandSide(yIn, dydx) ; // 1st Stage - Not doing, getting passed
274
275 for(i=0; i<numberOfVariables; ++i)
276 {
277 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
278 }
279 RightHandSide(yTemp, ak2) ; // 2nd Stage
280
281 for(i=0; i<numberOfVariables; ++i)
282 {
283 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
284 }
285 RightHandSide(yTemp, ak3) ; // 3rd Stage
286
287 for(i=0; i<numberOfVariables; ++i)
288 {
289 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
290 }
291 RightHandSide(yTemp, ak4) ; // 4th Stage
292
293 for(i=0; i<numberOfVariables; ++i)
294 {
295 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
296 b54*ak4[i]) ;
297 }
298 RightHandSide(yTemp, ak5) ; // 5th Stage
299
300 for(i=0; i<numberOfVariables; ++i)
301 {
302 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
303 b64*ak4[i] + b65*ak5[i]) ;
304 }
305 RightHandSide(yTemp, ak6) ; // 6th Stage
306
307 for(i=0; i<numberOfVariables; ++i)
308 {
309 yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
310 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
311 }
312 RightHandSide(yTemp, ak7); // 7th Stage
313
314 for(i=0; i<numberOfVariables; ++i)
315 {
316 yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
317 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
318 b87*ak7[i]);
319 }
320 RightHandSide(yTemp, ak8); // 8th Stage
321
322 for(i=0; i<numberOfVariables; ++i)
323 {
324 yOut[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
325 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
326 b97*ak7[i] + b98*ak8[i] );
327 }
328 RightHandSide(yOut, ak9); // 9th Stage
329
330
331 for(i=0; i<numberOfVariables; ++i)
332 {
333 // Estimate error as difference between 5th and
334 // 6th order methods
335 //
336 yErr[i] = Step*( dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
337 + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
338 + dc9*ak9[i] ) ;
339
340 // Saving 'estimated' derivative at end-point
341 // nextDydx[i] = ak9[i];
342
343 // Store Input and Final values, for possible use in calculating chord
344 //
345 fLastInitialVector[i] = yIn[i] ;
346 fLastFinalVector[i] = yOut[i];
347 fLastDyDx[i] = dydx[i];
348 }
349
350 fLastStepLength = Step;
351
352 return ;
353}

◆ StepperType()

G4StepperType G4DormandPrinceRK56::StepperType ( ) const
inlineoverridevirtual

Returns the stepper type-ID, "kDormandPrinceRK56".

Reimplemented from G4MagIntegratorStepper.

Definition at line 100 of file G4DormandPrinceRK56.hh.

100{ return kDormandPrinceRK56; }
kDormandPrinceRK56
@ kDormandPrinceRK56
G4DormandPrinceRK56.
Definition G4FieldParameters.hh:106
The documentation for this class was generated from the following files: