Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4DormandPrinceRK56 Class Reference

#include <G4DormandPrinceRK56.hh>

+ Inheritance diagram for G4DormandPrinceRK56:

Public Member Functions

 G4DormandPrinceRK56 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4DormandPrinceRK56 ()
 
 G4DormandPrinceRK56 (const G4DormandPrinceRK56 &)=delete
 
G4DormandPrinceRK56 & operator= (const G4DormandPrinceRK56 &)=delete
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
G4double DistChord () const
 
G4int IntegratorOrder () const
 
void SetupInterpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step)
 
void Interpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
 
void SetupInterpolation ()
 
void SetupInterpolate (const G4double yInput[], const G4double dydx[], const G4double Step)
 
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
 
virtual void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
 
virtual G4double DistChord () const =0
 
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
 
virtual G4int IntegratorOrder () const =0
 
G4int IntegrationOrder ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
const G4EquationOfMotion * GetEquationOfMotion () const
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
G4bool IsFSAL () const
 

Additional Inherited Members

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

Detailed Description

Definition at line 40 of file G4DormandPrinceRK56.hh.

Constructor & Destructor Documentation

◆ G4DormandPrinceRK56() [1/2]

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

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 ( )

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

Member Function Documentation

◆ DistChord()

G4double G4DormandPrinceRK56::DistChord ( ) const
virtual

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
G4DormandPrinceRK56::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4DormandPrinceRK56.cc:119
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:72

◆ IntegratorOrder()

G4int G4DormandPrinceRK56::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 60 of file G4DormandPrinceRK56.hh.

60{ return 5; }

◆ Interpolate() [1/2]

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

Definition at line 86 of file G4DormandPrinceRK56.hh.

91 {
92 Interpolate_low( yInput, dydx, Step, yOut, tau);
93 }
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 96 of file G4DormandPrinceRK56.hh.

97 {
98 Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau );
99 }
G4DormandPrinceRK56::Interpolate
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition: G4DormandPrinceRK56.hh:86

◆ Interpolate_high()

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

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 
)

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

Definition at line 79 of file G4DormandPrinceRK56.hh.

82 {
83 SetupInterpolate_low( yInput, dydx, Step);
84 }
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 
)

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 
)

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 74 of file G4DormandPrinceRK56.hh.

75 {
76 SetupInterpolate( fLastInitialVector, fLastDyDx, fLastStepLength);
77 }
G4DormandPrinceRK56::SetupInterpolate
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition: G4DormandPrinceRK56.hh:79

◆ Stepper()

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

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}

Referenced by DistChord().

The documentation for this class was generated from the following files: