RADAU5Run Method

public void Run(
	int N,
	IFVPOL FCN,
	ref double X,
	ref double[] Y,
	int offset_y,
	double XEND,
	ref double H,
	ref double[] RTOL,
	int offset_rtol,
	ref double[] ATOL,
	int offset_atol,
	int ITOL,
	IJVPOL JAC,
	int IJAC,
	ref int MLJAC,
	ref int MUJAC,
	IBBAMPL MAS,
	int IMAS,
	int MLMAS,
	ref int MUMAS,
	ISOLOUTR SOLOUT,
	int IOUT,
	ref double[] WORK,
	int offset_work,
	int LWORK,
	ref int[] IWORK,
	int offset_iwork,
	int LIWORK,
	double[] RPAR,
	int offset_rpar,
	int[] IPAR,
	int offset_ipar,
	ref int IDID
)

Public Sub Run ( 
	N As Integer,
	FCN As IFVPOL,
	ByRef X As Double,
	ByRef Y As Double(),
	offset_y As Integer,
	XEND As Double,
	ByRef H As Double,
	ByRef RTOL As Double(),
	offset_rtol As Integer,
	ByRef ATOL As Double(),
	offset_atol As Integer,
	ITOL As Integer,
	JAC As IJVPOL,
	IJAC As Integer,
	ByRef MLJAC As Integer,
	ByRef MUJAC As Integer,
	MAS As IBBAMPL,
	IMAS As Integer,
	MLMAS As Integer,
	ByRef MUMAS As Integer,
	SOLOUT As ISOLOUTR,
	IOUT As Integer,
	ByRef WORK As Double(),
	offset_work As Integer,
	LWORK As Integer,
	ByRef IWORK As Integer(),
	offset_iwork As Integer,
	LIWORK As Integer,
	RPAR As Double(),
	offset_rpar As Integer,
	IPAR As Integer(),
	offset_ipar As Integer,
	ByRef IDID As Integer
)

Parameters

N  Int32

DIMENSION OF THE SYSTEM

FCN  IFVPOL

NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE VALUE OF F(X,Y): SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) DOUBLE PRECISION X,Y(N),F(N) F(1)=... ETC. RPAR, IPAR (SEE BELOW)

X  Double

INITIAL X-VALUE

Y  Double

 

offset_y  Int32

 

XEND  Double

FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)

H  Double

INITIAL STEP SIZE GUESS; FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).

RTOL  Double

 

offset_rtol  Int32

 

ATOL  Double

 

offset_atol  Int32

 

ITOL  Int32

SWITCH FOR RTOL AND ATOL: ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF Y(I) BELOW RTOL*ABS(Y(I))+ATOL ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW RTOL(I)*ABS(Y(I))+ATOL(I).

JAC  IJVPOL

NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY A DUMMY SUBROUTINE IN THE CASE IJAC=0). FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR) DOUBLE PRECISION X,Y(N),DFY(LDFY,N) DFY(1,1)= ... LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS FURNISHED BY THE CALLING PROGRAM. IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO BE FULL AND THE PARTIAL DERIVATIVES ARE STORED IN DFY AS DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) ELSE, THE JACOBIAN IS TAKEN AS BANDED AND THE PARTIAL DERIVATIVES ARE STORED DIAGONAL-WISE AS DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).

IJAC  Int32

SWITCH FOR THE COMPUTATION OF THE JACOBIAN: IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.

MLJAC  Int32

SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 0.LE.MLJAC.LT.N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN MATRIX (.GE. NUMBER OF NON-ZERO DIAGONALS BELOW THE MAIN DIAGONAL).

MUJAC  Int32

UPPER BANDWITH OF JACOBIAN MATRIX (.GE. NUMBER OF NON- ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). NEED NOT BE DEFINED IF MLJAC=N.

MAS  IBBAMPL

NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- MATRIX M. IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY MATRIX AND NEEDS NOT TO BE DEFINED; SUPPLY A DUMMY SUBROUTINE IN THIS CASE. IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR) DOUBLE PRECISION AM(LMAS,N) AM(1,1)= .... IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED AS FULL MATRIX LIKE AM(I,J) = M(I,J) ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED DIAGONAL-WISE AS AM(I-J+MUMAS+1,J) = M(I,J).

IMAS  Int32

GIVES INFORMATION ON THE MASS-MATRIX: IMAS=0: M IS SUPPOSED TO BE THE IDENTITY MATRIX, MAS IS NEVER CALLED. IMAS=1: MASS-MATRIX IS SUPPLIED.

MLMAS  Int32

SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 0.LE.MLMAS.LT.N: MLMAS IS THE LOWER BANDWITH OF THE MATRIX (.GE. NUMBER OF NON-ZERO DIAGONALS BELOW THE MAIN DIAGONAL). MLMAS IS SUPPOSED TO BE .LE. MLJAC.

MUMAS  Int32

UPPER BANDWITH OF MASS-MATRIX (.GE. NUMBER OF NON- ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). NEED NOT BE DEFINED IF MLMAS=N. MUMAS IS SUPPOSED TO BE .LE. MUJAC.

SOLOUT  ISOLOUTR

NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE NUMERICAL SOLUTION DURING INTEGRATION. IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. SUPPLY A DUMMY SUBROUTINE IF IOUT=0. IT MUST HAVE THE FORM SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N, RPAR,IPAR,IRTRN) DOUBLE PRECISION X,Y(N),CONT(LRC) .... SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH GRID-POINT "X" (THEREBY THE INITIAL VALUE IS THE FIRST GRID-POINT). "XOLD" IS THE PRECEEDING GRID-POINT. "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN IS SET .LT.0, RADAU5 RETURNS TO THE CALLING PROGRAM. ----- CONTINUOUS OUTPUT: ----- DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH THE FUNCTION .GT..GT..GT. CONTR5(I,S,CONT,LRC) .LT..LT..LT. WHICH PROVIDES AN APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE S SHOULD LIE IN THE INTERVAL [XOLD,X]. DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE DENSE OUTPUT FUNCTION IS USED.

IOUT  Int32

SWITCH FOR CALLING THE SUBROUTINE SOLOUT: IOUT=0: SUBROUTINE IS NEVER CALLED IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.

WORK  Double

ARRAY OF WORKING SPACE OF LENGTH "LWORK". WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS FOR THE CODE. FOR STANDARD USE OF THE CODE WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE. WORK(21),..,WORK(LWORK) SERVE AS WORKING SPACE FOR ALL VECTORS AND MATRICES. "LWORK" MUST BE AT LEAST N*(LJAC+LMAS+3*LE+12)+20 WHERE LJAC=N IF MLJAC=N (FULL JACOBIAN) LJAC=MLJAC+MUJAC+1 IF MLJAC.LT.N (BANDED JAC.) AND LMAS=0 IF IMAS=0 LMAS=N IF IMAS=1 AND MLMAS=N (FULL) LMAS=MLMAS+MUMAS+1 IF MLMAS.LT.N (BANDED MASS-M.) AND LE=N IF MLJAC=N (FULL JACOBIAN) LE=2*MLJAC+MUJAC+1 IF MLJAC.LT.N (BANDED JAC.) IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM STORAGE REQUIREMENT IS LWORK = 4*N*N+12*N+20. IF IWORK(9)=M1.GT.0 THEN "LWORK" MUST BE AT LEAST N*(LJAC+12)+(N-M1)*(LMAS+3*LE)+20 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE NUMBER N CAN BE REPLACED BY N-M1.

offset_work  Int32

 

LWORK  Int32

DECLARED LENGTH OF ARRAY "WORK".

IWORK  Int32

INTEGER WORKING SPACE OF LENGTH "LIWORK". IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS FOR THE CODE. FOR STANDARD USE, SET IWORK(1),.., IWORK(20) TO ZERO BEFORE CALLING. IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA. "LIWORK" MUST BE AT LEAST 3*N+20.

offset_iwork  Int32

 

LIWORK  Int32

DECLARED LENGTH OF ARRAY "IWORK".

RPAR  Double

 

offset_rpar  Int32

 

IPAR  Int32

 

offset_ipar  Int32

 

IDID  Int32

REPORTS ON SUCCESSFULNESS UPON RETURN: IDID= 1 COMPUTATION SUCCESSFUL, IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT) IDID=-1 INPUT IS NOT CONSISTENT, IDID=-2 LARGER NMAX IS NEEDED, IDID=-3 STEP SIZE BECOMES TOO SMALL, IDID=-4 MATRIX IS REPEATEDLY SINGULAR.

Reference