Executes the Radau5 implicit Runge-Kutta solver for a stiff ODE system.
Namespace: Altaxo.Calc.Ode.Obsolete.Radau5Assembly: AltaxoCore (in AltaxoCore.dll) Version: 4.8.3572.0 (4.8.3572.0)
Syntaxpublic 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
)
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
- Initial values on entry and the computed solution vector on return.
- offset_y Int32
- The starting offset in Y.
- 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
- Relative tolerance array or scalar storage.
- offset_rtol Int32
- The starting offset in RTOL.
- ATOL Double
- Absolute tolerance array or scalar storage.
- offset_atol Int32
- The starting offset in ATOL.
- 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
- The starting offset in WORK.
- 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
- The starting offset in IWORK.
- LIWORK Int32
-
DECLARED LENGTH OF ARRAY "IWORK".
- RPAR Double
- User-defined real parameter array passed to callback routines.
- offset_rpar Int32
- The starting offset in RPAR.
- IPAR Int32
- User-defined integer parameter array passed to callback routines.
- offset_ipar Int32
- The starting offset in IPAR.
- 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.
See Also