/*************************************************************** * NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY * * DIFFERENTIAL EQUATIONS Y'=F(X,Y) BY ROSENBROCK METHOD. * * ------------------------------------------------------------ * * SAMPLE RUN: * * Example #1: * * (Solve set of DIFFERENTIAL equations (N=2): * * F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X) * * F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X)) * * Find values of F(1), F(2) at X=1.5]. * * * * SOLUTION AT X = 1.500000 * * Y(1) = 1.236003 * * Y(2) = -0.104930 * * * * LAST STEP SIZE = 0.000284 * * ERROR CODE = 1 * * * * Example #2: * * (Solve set of DIFFERENTIAL equations (N=5): * * F(1) = Y(2) * * F(2) = Y(3) * * F(3) = Y(4) * * F(4) = Y(5) * * F(5) = (45.0 * Y(3) * Y(4) * Y(5) - * * 40.0 * Y(4) * Y(4) * Y(4)) / (NINE * Y(3) * Y(3)) * * Find values of F(1), F(2), ..., F(5) at X=1.5). * * * * SOLUTION AT X = 1.500000 * * Y(1) = 4.363961 * * Y(2) = 4.000000 * * Y(3) = 2.828427 * * Y(4) = 0.000000 * * Y(5) = -3.771236 * * * * LAST STEP SIZE = 0.012570 * * ERROR CODE = 1 * * ------------------------------------------------------------ * * Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.0 By J-P Moreau, Paris * * (www.jpmoreau.fr) * **************************************************************** LIST OF USED FUNCTIONS (HERE INCLUDED): ======================================= ROS4, RO4COR, SHAMP, GRK4A, GRK4T, VELDD, VELDS, LSTAB, DECA, DECB, SOL, SOLB. Also uses files basis.h, basis_r.cpp, vmblock.h, vmblock.cpp. ---------------------------------------------------------------- Note: the banded matrix branch has not been tested here, but is fully implemented. */ #include #include #define NMX 25 typedef double MAT5[6][6]; //index 0 not used here double X,XEND, H,RTOL,ATOL; int I,IDID,N; double *Y, *WORK; int *IWORK; // variables for statistics (optional use) long NFCN,NSTEP,NJAC,NACCPT,NREJCT,NDEC,NSOL; int IDFX,IFCN,IJAC,IMAS,IOUT,ITOL,MLJAC,MLMAS; int LE1,LIWORK,LJAC,LMAS,LWORK,MUJAC,MUMAS; void *vmblock = NULL; /* define example #1 void FCN(int N, double X, double *Y, double *F) { F[1] = Y[1] * Y[2] + COS(X) - HALF * SIN(TWO * X); F[2] = Y[1] * Y[1] + Y[2] * Y[2] - (ONE + SIN(X)); } */ // define example #2 void FCN(int N, double X, double *Y, double *F) { F[1] = Y[2]; F[2] = Y[3]; F[3] = Y[4]; F[4] = Y[5]; F[5] = (45.0 * Y[3] * Y[4] * Y[5] - 40.0 * Y[4] * Y[4] * Y[4]) / (9.0 * Y[3] * Y[3]); } //Emulation of some Fortran intrinsic functions int IMAX(int a, int b) { if (a > b) return a; else return b; } int IMIN(int a, int b) { if (a < b) return a; else return b; } double XMAX(double a, double b) { if (a > b) return a; else return b; } double XMIN(double a, double b) { if (a < b) return a; else return b; } double XSIGN(double a, double b) { if (b<0) return (-ABS(a)); else return (ABS(a)); } //Header of core function called by ROS4. void RO4COR(int N, double *X, double *Y, double *XEND, double HMAX, double *H, double RTOL, double ATOL, int ITOL, int IJAC, int MLJAC, int MUJAC, int IDFX, int MLMAS, int MUMAS, int IOUT, int *IDID, long NMAX, double UROUND, int METH, double FAC1, double FAC2, double FACREJ, bool AUTNMS, bool IMPLCT, bool BANDED, int LFJAC, int LE, int LDMAS, double *YNEW, double *DY1, double *DY, double *AK1, double *AK2, double *AK3, double *AK4, double *FX, double *FJAC1, double *EE1, double *FMAS1, int *IP); //********************************************************************* void ROS4(int N,int IFCN,double *X, double *Y, double *XEND, double *H, double RTOL, double ATOL, int ITOL, int IJAC, int MLJAC, int MUJAC, int IDFX, int IMAS, int MLMAS, int MUMAS, int IOUT, double *WORK, int LWORK, int *IWORK, int LIWORK, int *IDID) { /* -------------------------------------------------------------------- { NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) { SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(X,Y). { THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 { WITH STEP SIZE CONTROL. { C.F. SECTION IV.7 { { AUTHORS: E. HAIRER AND G. WANNER { UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES { CH-1211 GENEVE 24, SWITZERLAND { E-MAIL: HAIRER@CGEUGE51.BITNET, WANNER@CGEUGE51.BITNET { { THIS CODE IS PART OF THE BOOK: { E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL { EQUATIONS II. STifF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. { SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, { SPRINGER-VERLAG [1990] { { VERSION OF OCTOBER 12, 1990 { { INPUT PARAMETERS { ---------------- { N DIMENSION OF THE SYSTEM { { FCN NAME OF PROCEDURE COMPUTING THE { VALUE OF F(X,Y): { Procedure FCN(N,X,Y,F) { X,Y[N],F[N]: Double; { F[1]=... ETC. { (Here removed from parameter list} { { IFCN GIVES INFORMATION ON FCN: { IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS) { IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS) { { X INITIAL X-VALUE { { Y[N] INITIAL VALUES FOR Y { { XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) { { H INITIAL STEP SIZE GUESS; { FOR STifF EQUATIONS WITH INITIAL TRANSIENT, { H=1.D0/(NORM OF F'), USUALLY 1e-2 OR 1e-3, IS GOOD. { THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY { ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW { STEPS IN PROCEDURE SOLOUT, WHEN YOU ARE NOT SURE. { (if H=0.0, THE CODE PUTS H=1e-6). { { RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY { CAN BE BOTH SCALARS OR else BOTH VECTORS OF LENGTH N. { (multiple tolerances not implemented here). { { ITOL 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 NAME OF THE PROCEDURE WHICH COMPUTES { THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y { (THIS ROUTINE IS ONLY CALLED if IJAC=1. { FOR IJAC=1, THIS PROCEDURE MUST HAVE THE FORM: { Procedure JAC(N,X,Y,DFY,LDFY) { X,Y[N],DFY[LDFY,N]:Double; { DFY[1,1]= ... { LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS { FURNISHED BY THE CALLING PROGRAM. { if MLJAC = 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 SWITCH FOR THE COMPUTATION OF THE JACOBIAN: { IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE { DifFERENCES, PROCEDURE JAC IS NEVER CALLED. { IJAC=1: JACOBIAN IS SUPPLIED BY PROCEDURE JAC. { { MLJAC 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<=MLJAC= NUMBER OF NON-ZERO DIAGONALS BELOW { THE MAIN DIAGONAL). { { MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- { ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). { NEED NOT BE DEFINED if MLJAC=N. { { DFX NAME OF THE PROCEDURE WHICH COMPUTES { THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X { THIS ROUTINE IS ONLY CALLED if IDFX=1 AND IFCN=1. { OTHERWISE, THIS PROCEDURE MUST HAVE THE FORM: { Procedure DFX(N,X,Y,FX) { X,Y[N],FX[N]:Double; { FX[1]= ... { { IDFX SWITCH FOR THE COMPUTATION OF THE DF/DX: { IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE { DifFERENCES, PROCEDURE DFX IS NEVER CALLED. { IDFX=1: DF/DX IS SUPPLIED BY PROCEDURE DFX. { { ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- { ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - { { MAS NAME OF PROCEDURE COMPUTING THE MASS-MATRIX M. { if IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY { MATRIX AND NEEDS NOT TO BE DEFINED. { if IMAS=1, THE PROCEDURE MAS IS OF THE FORM: { Procedure MAS(N,AM,LMAS) { AM(LMAS,N):Double; { AM[1,1]= .... { if MLMAS = 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 GIVES INFORMATION ON THE MASS-MATRIX: { IMAS=0: M IS SUPPOSED TO BE THE IDENTITY { MATRIX, PROCEDURE MAS IS NEVER CALLED. { IMAS=1: MASS-MATRIX IS SUPPLIED. { { MLMAS 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<=MLMAS= NUMBER OF NON-ZERO DIAGONALS BELOW { THE MAIN DIAGONAL). { MLMAS IS SUPPOSED TO BE <= MLJAC. { { MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- { ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). { NEED NOT BE DEFINED if MLMAS=N. { MUMAS IS SUPPOSED TO BE <= MUJAC. { { SOLOUT NAME OF PROCEDURE PROVIDING THE { NUMERICAL SOLUTION DURING INTEGRATION. { if IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. { IT MUST HAVE THE FORM: { Procedure SOLOUT (NR,XOLD,X,Y,N,IRTRN) { X,Y[N]:Double; { .... { SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH { GRID-POINT "X" (THEREBY THE INITIAL VALUE IS { THE FIRST GRID-POINT). { "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. if IRTRN { IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM. { { IOUT GIVES INFORMATION ON THE PROCEDURE SOLOUT: { IOUT=0: Procedure IS NEVER CALLED { IOUT=1: Procedure IS USED FOR OUTPUT { { WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". { SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES. { "LWORK" MUST BE AT LEAST { N*[LJAC+LMAS+LE1+8]+5 { WHERE { LJAC=N if MLJAC=N (FULL JACOBIAN) { LJAC=MLJAC+MUJAC+1 if MLJAC= 7) { printf(" CURIOUS INPUT IWORK[2]=%d\n", IWORK[2]); ARRET=TRUE; } } // -------- UROUND SMALLEST NUMBER SATISFYING 1 + UROUND > 1 ----- if (WORK[1] == ZERO) UROUND=1e-16; else { UROUND=WORK[1]; if (UROUND <= 1e-14 || UROUND >= ONE) { printf(" COEFFICIENTS HAVE 16 DIGITS, UROUND=%f\n", WORK[1]); ARRET=TRUE; } } // -------- MAXIMAL STEP SIZE -------- if (WORK[2] == ZERO) HMAX=XEND-X; else HMAX=WORK[2]; // ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION ------- if (WORK[3] == ZERO) FAC1=5.0; else FAC1=ONE/WORK[3]; if (WORK[4] == ZERO) FAC2=ONE/6.0; else FAC2=ONE/WORK[4]; // ------- FACREJ FOR THE HUMP ------ if (WORK[5] == ZERO) FACREJ=0.1; else FACREJ=WORK[5]; // --------- CHECK if TOLERANCES ARE O.K. ------ if (ITOL == 0) if (ATOL <= ZERO || RTOL <= 10.0*UROUND) { printf(" TOLERANCES ARE TOO SMALL\n"); ARRET=TRUE; } else { /* Multiple tolerances not implemented here for (I=1; I<=N; I++) if (ATOL[I] <= ZERO || RTOL[I] <= 10.0*UROUND) { printf(" TOLERANCES(%d) ARE TOO SMALL\n", I); ARRET=TRUE; } */ } /*** *** *** *** *** *** *** *** *** *** *** *** *** COMPUTATION OF ARRAY ENTRIES *** *** *** *** *** *** *** *** *** *** *** *** *** ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ? */ AUTNMS=(IFCN==0); IMPLCT=(IMAS!=0); JBAND=(MLJAC!=N); ARRET=FALSE; // --- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- // --- JACOBIAN --- if (JBAND) LDJAC=MLJAC+MUJAC+1; else LDJAC=N; // -- MATRIX E FOR LINEAR ALGEBRA -- if (JBAND) LDE=2*MLJAC+MUJAC+1; else LDE=N; // -- MASS MATRIX -- if (IMPLCT) { if (MLMAS != N) LDMAS=MLMAS+MUMAS+1; else LDMAS=N; // ------ BANDWITH OF MAS MUST NOT BE LARGER THAN BANDWITH OF JAC ---- if (MLMAS > MLJAC || MUMAS > MUJAC) { printf(" BANDWITH OF MAS NOT LARGER THAN BANDWITH OF JAC\n"); ARRET=TRUE; } } else LDMAS=0; LDMAS2=IMAX(1,LDMAS); // ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- IEYNEW=6; IEDY1=IEYNEW+N; IEDY=IEDY1+N; IEAK1=IEDY+N; IEAK2=IEAK1+N; IEAK3=IEAK2+N; IEAK4=IEAK3+N; IEFX =IEAK4+N; IEJAC=IEFX +N; IEMAS=IEJAC+N*LDJAC; IEE =IEMAS+N*LDMAS; // -------- TOTAL STORAGE REQUIREMENT ----------- ISTORE=IEE+N*LDE-1; if (ISTORE > LWORK) { printf(" INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=%d\n", ISTORE); ARRET=TRUE; } // ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- IEIP=3; // ------------ TOTAL REQUIREMENT --------------- ISTORE=IEIP+N-1; if (ISTORE > LIWORK) { printf(" INSUFF. STORAGE FOR IWORK, MIN. LIWORK=%d\n", ISTORE); ARRET=TRUE; } // ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 ----- if (ARRET) { *IDID=-1; goto e10; } // --- prepare arguments of RO4COR --- for (I=1; I<=N; I++) if (I+IEYNEW-1 <= N) TMP1[I] = WORK[I+IEYNEW-1]; else TMP1[I] = ZERO; for (I=1; I<=N; I++) if (I+IEDY1-1 <= N) TMP2[I] = WORK[I+IEDY1-1]; else TMP2[I] = ZERO; for (I=1; I<=N; I++) if (I+IEDY-1 <= N) TMP3[I] = WORK[I+IEDY-1]; else TMP3[I] = ZERO; for (I=1; I<=N; I++) if (I+IEAK1-1 <= N) TMP4[I] = WORK[I+IEAK1-1]; else TMP4[I] = ZERO; for (I=1; I<=N; I++) if (I+IEAK2-1 <= N) TMP5[I] = WORK[I+IEAK2-1]; else TMP5[I] = ZERO; for (I=1; I<=N; I++) if (I+IEAK3-1 <= N) TMP6[I] = WORK[I+IEAK3-1]; else TMP6[I] = ZERO; for (I=1; I<=N; I++) if (I+IEAK4-1 <= N) TMP7[I] = WORK[I+IEAK4-1]; else TMP7[I] = ZERO; for (I=1; I<=N; I++) if (I+IEFX-1 <= N) TMP8[I] = WORK[I+IEFX-1]; else TMP8[I] = ZERO; for (I=1; I<=N; I++) if (I+IEJAC-1 <= N) TMP9[I] = WORK[I+IEJAC-1]; else TMP9[I] = ZERO; for (I=1; I<=N; I++) if (I+IEE-1 <= N) TMP10[I] = WORK[I+IEE-1]; else TMP10[I] = ZERO; for (I=1; I<=N; I++) if (I+IEMAS-1 <= N) TMP11[I] = WORK[I+IEMAS-1]; else TMP11[I] = ZERO; for (I=1; I<=N; I++) if (I+IEIP-1 <= N) ITMP[I] = IWORK[I+IEIP-1]; else ITMP[I] = 0; // ---------- CALL TO CORE INTEGRATOR --------------- RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IJAC, MLJAC,MUJAC,IDFX,MLMAS,MUMAS,IOUT,IDID, NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS, IMPLCT,JBAND,LDJAC,LDE,LDMAS2,TMP1,TMP2, TMP3,TMP4,TMP5,TMP6,TMP7,TMP8,TMP9,TMP10, TMP11, ITMP); // ----------- RETURN SECTION ----------- e10: vmfree(vmblock); } //ROS4() //Headers of functions used below void SHAMP(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void GRK4T(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void GRK4A(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void VELDD(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void VELDS(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void LSTAB(double *A21,double *A31,double *A32,double *C21,double *C31, double *C32,double *C41,double *C42,double *C43,double *B1, double *B2,double *B3,double *B4,double *E1,double *E2, double *E3,double *E4,double *GAMMA,double *C2, double *C3,double *D1,double *D2,double *D3,double *D4); void DECA(int, int, MAT5, int *, int *); void SOL (int, int, MAT5, double *, int *); void DECB(int, int, MAT5, int, int, int *, int *); void SOLB(int, int, MAT5, int, int, double *, int *); // --------- ... AND HERE IS THE CORE INTEGRATOR ---------- void RO4COR(int N, double *X, double *Y, double *XEND, double HMAX, double *H, double RTOL, double ATOL,int ITOL,int IJAC, int MLJAC,int MUJAC, int IDFX,int MLMAS,int MUMAS,int IOUT, int *IDID, long NMAX, double UROUND, int METH, double FAC1, double FAC2,double FACREJ, bool AUTNMS,bool IMPLCT,bool BANDED, int LFJAC,int LE,int LDMAS, double *YNEW,double *DY1,double *DY, double *AK1,double *AK2,double *AK3, double *AK4,double *FX, double *FJAC1,double *EE1,double *FMAS1, int *IP) { /* --------------------------------------------------------- { CORE INTEGRATOR FOR ROS4 { PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED { ---------------------------------------------------------- { DECLARATIONS { ---------------------------------------------------------*/ // Labels: e1,e2, e12, e14, e79, e80 int I,J,K,L,MBB,MBDIAG,MBJAC,MDIAG,MDIFF,MLE,MUE; bool REJECT,RJECT2; MAT5 FJAC, E, FMAS; double A21,A31,A32,C21,C31,C32,C41,C42,C43; double B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4; double DELT,HMAXN,HOPT,POSNEG,XDELT,XXOLD,YSAFE; int IRTRN,LBEG,LEND,MD,MUJACJ,MUJACP,NSING; double FAC,HC21,HC31,HC32,HC41,HC42,HC43; int I1,I2,INFO; double ERR, HD1,HD2,HD3,HD4, HNEW,S,SK; // ------- restore parameters FJAC, E, FMAS ---------- for (J=1; J<=N; J++) for (I=1; I<=N; I++) { FJAC[I][J] = FJAC1[I+J-1]; E[I][J] = EE1[I+J-1]; FMAS[I][J] = FMAS1[I+J-1]; } /* ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- if (IMPLCT) MAS(N,FMAS,LDMAS); ---- PREPARE BANDWIDTHS ----- */ if (BANDED) { MLE=MLJAC; MUE=MUJAC; MBJAC=MLJAC+MUJAC+1; MBB=MLMAS+MUMAS+1; MDIAG=MLE+MUE+1; MBDIAG=MUMAS+1; MDIFF=MLE+MUE-MUMAS; } /*** *** *** *** *** *** *** INITIALISATIONS *** *** *** *** *** *** ***/ POSNEG=XSIGN(ONE,*XEND-*X); //Bug corrected 01/05/2012 if (METH == 1) SHAMP(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); if (METH == 2) GRK4T(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); if (METH == 3) GRK4A(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); if (METH == 4) VELDS(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); if (METH == 5) VELDD(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); if (METH == 6) LSTAB(&A21,&A31,&A32,&C21,&C31,&C32,&C41,&C42,&C43, &B1,&B2,&B3,&B4,&E1,&E2,&E3,&E4,&GAMMA, &C2,&C3,&D1,&D2,&D3,&D4); // --- INITIAL PREPARATIONS --- HMAXN=XMIN(ABS(HMAX),ABS(*XEND-*X)); *H=XMIN(XMAX(1e-10,ABS(*H)),HMAXN); *H=XSIGN(*H,POSNEG); REJECT=FALSE; NSING=0; IRTRN=1; XXOLD=*X; if (IOUT != 0) //SOLOUT(NACCPT+1,XXOLD,X,Y,N,IRTRN); if (IRTRN < 0) goto e79; // --- BASIC INTEGRATION STEP --- e1: if (NSTEP > NMAX || *X+0.1*(*H) == *X ||ABS(*H) <= UROUND) goto e79; if ((*X-*XEND)*POSNEG+UROUND > ZERO) { *H=HOPT; *IDID=1; return; } HOPT=*H; if ((*X+*H-*XEND)*POSNEG > ZERO) *H=(*XEND-*X); FCN(N,*X,Y,DY1); NFCN++; /*** *** *** *** *** *** *** COMPUTATION OF THE JACOBIAN *** *** *** *** *** *** ***/ NJAC++; if (IJAC == 0) { // --- COMPUTE JACOBIAN MATRIX NUMERICALLY --- if (BANDED) { // --- JACOBIAN IS BANDED --- MUJACP=MUJAC+1; MD=IMIN(MBJAC,N); for (K=1; K<=MD; K++) { J=K; e12: AK2[J]=Y[J]; AK3[J]=SQRT(UROUND*XMAX(1e-5,ABS(Y[J]))); Y[J]=Y[J]+AK3[J]; J=J+MD; if (J <= N) goto e12; FCN(N,*X,Y,AK1); J=K; LBEG=IMAX(1,J-MUJAC); e14: LEND=IMIN(N,J+MLJAC); Y[J]=AK2[J]; MUJACJ=MUJACP-J; for (L=LBEG; L<=LEND; L++) FJAC[L+MUJACJ][J]=(AK1[L]-DY1[L])/AK3[J]; J=J+MD; LBEG=LEND+1; if (J <= N) goto e14; } } else { // --- JACOBIAN IS FULL --- for (I=1; I<=N; I++) { YSAFE=Y[I]; DELT=SQRT(UROUND*XMAX(1e-5,ABS(YSAFE))); Y[I]=YSAFE+DELT; FCN(N,*X,Y,AK1); for (J=1; J<=N; J++) FJAC[J][I]=(AK1[J]-DY1[J])/DELT; Y[I]=YSAFE; } MLJAC=N; } } else // --- COMPUTE JACOBIAN MATRIX ANALYTICALLY --- // JAC(N,X,Y,FJAC,LFJAC); if (!AUTNMS) if (IDFX == 0) { // --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X --- DELT=SQRT(UROUND*XMAX(1e-5,ABS(*X))); XDELT=*X+DELT; FCN(N,XDELT,Y,AK1); for (J=1; J<=N; J++) FX[J]=(AK1[J]-DY1[J])/DELT; } else // --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X --- // DFX(N,X,Y,FX); /*** *** *** *** *** *** *** COMPUTE THE STAGES *** *** *** *** *** *** ***/ e2: NDEC++; HC21=C21/(*H); HC31=C31/(*H); HC32=C32/(*H); HC41=C41/(*H); HC42=C42/(*H); HC43=C43/(*H); FAC=ONE/((*H)*GAMMA); // Note: This part corresponds to IMPLCT=0, (THE DIFFERENTIAL EQUATION // IS IN EXPLICIT FORM). The part for IMPLCT<>0 is not implemented here. if (BANDED) { // --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) for (J=1; J<=N; J++) { I1=IMAX(1,MUJAC+2-J); I2=IMIN(MBJAC,N+MUJAC+1-J); for (I=I1; I<=I2; I++) E[I+MLE][J]=-FJAC[I][J]; E[MDIAG][J] += FAC; } DECB(N,LE,E,MLE,MUE,IP,&INFO); if (INFO != 0) goto e80; if (AUTNMS) { /*--- THIS PART COMPUTES THE STAGES IN THE CASE WHERE --- 1] THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM --- 2] THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX --- 3] THE DIFFERENTIAL EQUATION IS AUTONOMOUS */ for (I=1; I<=N; I++) AK1[I]=DY1[I]; SOLB(N,LE,E,MLE,MUE,AK1,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A21*AK1[I]; FCN(N,*X,YNEW,DY); for (I=1; I<=N; I++) AK2[I]=DY[I]+HC21*AK1[I]; SOLB(N,LE,E,MLE,MUE,AK2,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A31*AK1[I]+A32*AK2[I]; FCN(N,*X,YNEW,DY); for (I=1; I<=N; I++) AK3[I]=DY[I]+HC31*AK1[I]+HC32*AK2[I]; SOLB(N,LE,E,MLE,MUE,AK3,IP); for (I=1; I<=N; I++) AK4[I]=DY[I]+HC41*AK1[I]+HC42*AK2[I]+HC43*AK3[I]; SOLB(N,LE,E,MLE,MUE,AK4,IP); } else { /*--- THIS PART COMPUTES THE STAGES IN THE CASE WHERE --- 1] THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM --- 2] THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX --- 3] THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS */ HD1=(*H)*D1; HD2=(*H)*D2; HD3=(*H)*D3; HD4=(*H)*D4; for (I=1; I<=N; I++) AK1[I]=DY1[I]+HD1*FX[I]; SOLB(N,LE,E,MLE,MUE,AK1,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A21*AK1[I]; FCN(N,*X+C2*(*H),YNEW,DY); for (I=1; I<=N; I++) AK2[I]=DY[I]+HD2*FX[I]+HC21*AK1[I]; SOLB(N,LE,E,MLE,MUE,AK2,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A31*AK1[I]+A32*AK2[I]; FCN(N,*X+C3*(*H),YNEW,DY); for (I=1; I<=N; I++) AK3[I]=DY[I]+HD3*FX[I]+HC31*AK1[I]+HC32*AK2[I]; SOLB(N,LE,E,MLE,MUE,AK3,IP); for (I=1; I<=N; I++) AK4[I]=DY[I]+HD4*FX[I]+HC41*AK1[I]+HC42*AK2[I]+HC43*AK3[I]; SOLB(N,LE,E,MLE,MUE,AK4,IP); } } else { // --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) for (J=1; J<=N; J++) { for (I=1; I<=N; I++) E[I][J]=-FJAC[I][J]; E[J][J] += FAC; } DECA(N,LE,E,IP,&INFO); if (INFO != 0) goto e80; if (AUTNMS) { /*--- THIS PART COMPUTES THE STAGES IN THE CASE WHERE --- 1] THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM --- 2] THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX --- 3] THE DIFFERENTIAL EQUATION IS AUTONOMOUS */ for (I=1; I<=N; I++) AK1[I]=DY1[I]; SOL(N,LE,E,AK1,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A21*AK1[I]; FCN(N,*X,YNEW,DY); for (I=1; I<=N; I++) AK2[I]=DY[I]+HC21*AK1[I]; SOL(N,LE,E,AK2,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A31*AK1[I]+A32*AK2[I]; FCN(N,*X,YNEW,DY); for (I=1; I<=N; I++) AK3[I]=DY[I]+HC31*AK1[I]+HC32*AK2[I]; SOL(N,LE,E,AK3,IP); for (I=1; I<=N; I++) AK4[I]=DY[I]+HC41*AK1[I]+HC42*AK2[I]+HC43*AK3[I]; SOL(N,LE,E,AK4,IP); } else { /*--- THIS PART COMPUTES THE STAGES IN THE CASE WHERE --- 1] THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM --- 2] THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX --- 3] THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS */ HD1=(*H)*D1; HD2=(*H)*D2; HD3=(*H)*D3; HD4=(*H)*D4; for (I=1; I<=N; I++) AK1[I]=DY1[I]+HD1*FX[I]; SOL(N,LE,E,AK1,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A21*AK1[I]; FCN(N,*X+C2*(*H),YNEW,DY); for (I=1; I<=N; I++) AK2[I]=DY[I]+HD2*FX[I]+HC21*AK1[I]; SOL(N,LE,E,AK2,IP); for (I=1; I<=N; I++) YNEW[I]=Y[I]+A31*AK1[I]+A32*AK2[I]; FCN(N,*X+C3*(*H),YNEW,DY); for (I=1; I<=N; I++) AK3[I]=DY[I]+HD3*FX[I]+HC31*AK1[I]+HC32*AK2[I]; SOL(N,LE,E,AK3,IP); for (I=1; I<=N; I++) AK4[I]=DY[I]+HD4*FX[I]+HC41*AK1[I]+HC42*AK2[I]+HC43*AK3[I]; SOL(N,LE,E,AK4,IP); } } NSOL=NSOL+4; NFCN=NFCN+2; /*** *** *** *** *** *** *** ERROR ESTIMATION *** *** *** *** *** *** ***/ NSTEP++; // ------------ NEW SOLUTION --------------- for (I=1; I<=N; I++) YNEW[I]=Y[I]+B1*AK1[I]+B2*AK2[I]+B3*AK3[I]+B4*AK4[I]; // ------------ COMPUTE ERROR ESTIMATION ---------------- ERR=ZERO; for (I=1; I<=N; I++) { S=E1*AK1[I]+E2*AK2[I]+E3*AK3[I]+E4*AK4[I]; if (ITOL == 0) SK=ATOL+RTOL*XMAX(ABS(Y[I]),ABS(YNEW[I])); else //SK=ATOL[I]+RTOL[I]*XMAX1(ABS(Y[I]),ABS(YNEW[I])); printf(" S=%f SK=%f\n", S, SK); ERR += (S/SK)*(S/SK); } ERR=SQRT(ERR/(1.0*N)); // --- COMPUTATION OF HNEW // --- WE REQUIRE 0.2<=HNEW/H<=6.0 FAC=XMAX(FAC2,XMIN(FAC1,pow(ERR,(0.25/0.9)))); HNEW=*H/FAC; /*** *** *** *** *** *** *** IS THE ERROR SMALL ENOUGH ? *** *** *** *** *** *** ***/ if (ERR <= ONE) { // --- STEP IS ACCEPTED --- NACCPT++; for (I=1; I<=N; I++) Y[I]=YNEW[I]; XXOLD=*X; *X=*X+*H; if (IOUT != 0) //SOLOUT(NACCPT+1,XXOLD,X,Y,N,IRTRN) if (IRTRN < 0) goto e79; if (ABS(HNEW) > HMAXN) HNEW=POSNEG*HMAXN; if (REJECT) HNEW=POSNEG*XMIN(ABS(HNEW),ABS(*H)); REJECT=FALSE; RJECT2=FALSE; *H=HNEW; goto e1; } else { // --- STEP IS REJECTED --- if (RJECT2) HNEW=*H*FACREJ; if (REJECT) RJECT2=TRUE; REJECT=TRUE; *H=HNEW; if (NACCPT >= 1) NREJCT++; goto e2; } // --- EXIT SECTION --- e80: printf(" MATRIX E IS SINGULAR, INFO = %d\n",INFO); NSING++; if (NSING >= 5) goto e79; *H=(*H)*HALF; goto e2; e79: printf("\n EXIT OF ROS4 AT X=%f H=%f\n", X, H); *IDID=-1; } //RO4COR() void SHAMP(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { *A21=2.0; *A31=48.0/25.0; *A32=6.0/25.0; *C21=-8.0; *C31=372.0/25.0; *C32=12.0/5.0; *C41=-112.0/125.0; *C42=-54.0/125.0; *C43=-2.0/5.0; *B1=19.0/9.0; *B2=1.0/2.0; *B3=25.0/108.0; *B4=125.0/108.0; *E1=17.0/54.0; *E2=7.0/36.0; *E3=0.0; *E4=125.0/108.0; *GAMMA=0.5; *C2= 0.1000000000000000e+01; *C3= 0.6000000000000000e+00; *D1= 0.5000000000000000e+00; *D2=-0.1500000000000000e+01; *D3= 0.2420000000000000e+01; *D4= 0.1160000000000000e+00; } void GRK4A(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { *A21= 0.1108860759493671e+01; *A31= 0.2377085261983360e+01; *A32= 0.1850114988899692e+00; *C21=-0.4920188402397641e+01; *C31= 0.1055588686048583e+01; *C32= 0.3351817267668938e+01; *C41= 0.3846869007049313e+01; *C42= 0.3427109241268180e+01; *C43=-0.2162408848753263e+01; *B1= 0.1845683240405840e+01; *B2= 0.1369796894360503e+00; *B3= 0.7129097783291559e+00; *B4= 0.6329113924050632e+00; *E1= 0.4831870177201765e-01; *E2=-0.6471108651049505e+00; *E3= 0.2186876660500240e+00; *E4=-0.6329113924050632e+00; *GAMMA= 0.3950000000000000e+00; *C2= 0.4380000000000000e+00; *C3= 0.8700000000000000e+00; *D1= 0.3950000000000000e+00; *D2=-0.3726723954840920e+00; *D3= 0.6629196544571492e-01; *D4= 0.4340946962568634e+00; } void GRK4T(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { *A21= 0.2000000000000000e+01; *A31= 0.4524708207373116e+01; *A32= 0.4163528788597648e+01; *C21=-0.5071675338776316e+01; *C31= 0.6020152728650786e+01; *C32= 0.1597506846727117e+00; *C41=-0.1856343618686113e+01; *C42=-0.8505380858179826e+01; *C43=-0.2084075136023187e+01; *B1= 0.3957503746640777e+01; *B2= 0.4624892388363313e+01; *B3= 0.6174772638750108e+00; *B4= 0.1282612945269037e+01; *E1= 0.2302155402932996e+01; *E2= 0.3073634485392623e+01; *E3=-0.8732808018045032e+00; *E4=-0.1282612945269037e+01; *GAMMA= 0.2310000000000000e+00; *C2= 0.4620000000000000e+00; *C3= 0.8802083333333334e+00; *D1= 0.2310000000000000e+00; *D2=-0.3962966775244303e-01; *D3= 0.5507789395789127e+00; *D4=-0.5535098457052764e-01; } void VELDS(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { // --- METHOD GIVEN BY VAN VELDHUIZEN --- *A21= 0.2000000000000000e+01; *A31= 0.1750000000000000e+01; *A32= 0.2500000000000000e+00; *C21=-0.8000000000000000e+01; *C31=-0.8000000000000000e+01; *C32=-0.1000000000000000e+01; *C41= 0.5000000000000000e+00; *C42=-0.5000000000000000e+00; *C43= 0.2000000000000000e+01; *B1= 0.1333333333333333e+01; *B2= 0.6666666666666667e+00; *B3=-0.1333333333333333e+01; *B4= 0.1333333333333333e+01; *E1=-0.3333333333333333e+00; *E2=-0.3333333333333333e+00; *E3=-0.0000000000000000e+00; *E4=-0.1333333333333333e+01; *GAMMA= 0.5000000000000000e+00; *C2= 0.1000000000000000e+01; *C3= 0.5000000000000000e+00; *D1= 0.5000000000000000e+00; *D2=-0.1500000000000000e+01; *D3=-0.7500000000000000e+00; *D4= 0.2500000000000000e+00; } void VELDD(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { // --- METHOD GIVEN BY VAN VELDHUIZEN --- *A21= 0.2000000000000000e+01; *A31= 0.4812234362695436e+01; *A32= 0.4578146956747842e+01; *C21=-0.5333333333333331e+01; *C31= 0.6100529678848254e+01; *C32= 0.1804736797378427e+01; *C41=-0.2540515456634749e+01; *C42=-0.9443746328915205e+01; *C43=-0.1988471753215993e+01; *B1= 0.4289339254654537e+01; *B2= 0.5036098482851414e+01; *B3= 0.6085736420673917e+00; *B4= 0.1355958941201148e+01; *E1= 0.2175672787531755e+01; *E2= 0.2950911222575741e+01; *E3=-0.7859744544887430e+00; *E4=-0.1355958941201148e+01; *GAMMA= 0.2257081148225682e+00; *C2= 0.4514162296451364e+00; *C3= 0.8755928946018455e+00; *D1= 0.2257081148225682e+00; *D2=-0.4599403502680582e-01; *D3= 0.5177590504944076e+00; *D4=-0.3805623938054428e-01; } void LSTAB(double *A21,double *A31,double *A32,double *C21,double *C31,double *C32, double *C41,double *C42,double *C43,double *B1,double *B2,double *B3, double *B4,double *E1,double *E2,double *E3,double *E4,double *GAMMA, double *C2,double *C3,double *D1,double *D2,double *D3,double *D4) { // --- AN L-STABLE METHOD --- *A21= 0.2000000000000000e+01; *A31= 0.1867943637803922e+01; *A32= 0.2344449711399156e+00; *C21=-0.7137615036412310e+01; *C31= 0.2580708087951457e+01; *C32= 0.6515950076447975e+00; *C41=-0.2137148994382534e+01; *C42=-0.3214669691237626e+00; *C43=-0.6949742501781779e+00; *B1= 0.2255570073418735e+01; *B2= 0.2870493262186792e+00; *B3= 0.4353179431840180e+00; *B4= 0.1093502252409163e+01; *E1=-0.2815431932141155e+00; *E2=-0.7276199124938920e-01; *E3=-0.1082196201495311e+00; *E4=-0.1093502252409163e+01; *GAMMA= 0.5728200000000000e+00; *C2= 0.1145640000000000e+01; *C3= 0.6552168638155900e+00; *D1= 0.5728200000000000e+00; *D2=-0.1769193891319233e+01; *D3= 0.7592633437920482e+00; *D4=-0.1049021087100450e+00; } void DECA(int N, int NDIM, MAT5 A, int *IP, int *IER) { // VERSION REAL DOUBLE PRECISION OF DEC // Labels: e20, e50, e70, e80 int NM1,K,KP1,M,I,J; double T; /*---------------------------------------------------------------------- { MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. { INPUTS: { N = ORDER OF MATRIX. { NDIM = DECLARED DIMENSION OF ARRAY A . { A = MATRIX TO BE TRIANGULARIZED. { OUTPUTS: { A[I,J], I <= J = UPPER TRIANGULAR FACTOR, U. { A[I,J], I > J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. { IP[K], K < N = INDEX OF K-TH PIVOT ROW. { IP[N] = (-1)^(NUMBER OF INTERCHANGES) OR O (^ for Power). { IER = 0 if MATRIX A IS NONSINGULAR, OR K if FOUND TO BE { SINGULAR AT STAGE K. { USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. { DETERM[A] = IP[N]*A[1,1]*A[2,2]*...*A[N,N]. { if IP[N]=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. { { REFERENCE.. { C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, { C.A.C.M. 15 [1972], P. 274. {----------------------------------------------------------------------*/ *IER = 0; IP[N] = 1; if (N == 1) goto e70; NM1 = N - 1; for (K=1; K<=NM1; K++) { KP1 = K + 1; M = K; for (I=KP1; I<=N; I++) if (ABS(A[I][K]) > ABS(A[M][K])) M = I; IP[K] = M; T = A[M][K]; if (M == K) goto e20; IP[N] = -IP[N]; A[M][K] = A[K][K]; A[K][K] = T; e20: if (T == ZERO) goto e80; T = ONE / T; for (I=KP1; I<=N; I++) A[I][K] = -A[I][K]*T; for (J=KP1; J<=N; J++) { T = A[M][J]; A[M][J] = A[K][J]; A[K][J] = T; if (T == ZERO) goto e50; for (I=KP1; I<=N; I++) A[I][J] += A[I][K]*T; e50: ;} } e70: K = N; if (A[N][N] == ZERO) goto e80; return; e80: *IER = K; IP[N] = 0; } // DECA() void SOL (int N, int NDIM, MAT5 A, double *B, int *IP) { // VERSION REAL DOUBLE PRECISION // Label: e50; int NM1,K,KP1,M,I,KB,KM1; double T; /*---------------------------------------------------------------------- { SOLUTION OF LINEAR SYSTEM, A*X = B. { INPUTS: { N = ORDER OF MATRIX. { NDIM = DECLARED DIMENSION OF ARRAY A. { A = TRIANGULARIZED MATRIX OBTAINED FROM DEC. { B = RIGHT HAND SIDE VECTOR. { IP = PIVOT VECTOR OBTAINED FROM DECA. { DO NOT USE if DECA HAS SET IER <> 0. { OUTPUT: { B = SOLUTION VECTOR, X. {----------------------------------------------------------------------*/ if (N == 1) goto e50; NM1 = N - 1; for (K=1; K<=NM1; K++) { KP1 = K + 1; M = IP[K]; T = B[M]; B[M] = B[K]; B[K] = T; for (I=KP1; I<=N; I++) B[I] += A[I][K]*T; } for (KB=1; KB<=NM1; KB++) { KM1 = N - KB; K = KM1 + 1; B[K] /= A[K][K]; T = -B[K]; for (I=1; I<=KM1; I++) B[I] += A[I][K]*T; } e50: B[1] /= A[1][1]; } void DECB(int N, int NDIM, MAT5 A, int ML, int MU, int *IP, int *IER) { // Labels: e7,e20,e30,e45, e60,e70,e80 double T; int I,IJK,J,JK,JU,K,KP1,M,MD,MD1,MDL,MM,NM1; /*---------------------------------------------------------------------- { MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED { MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU { INPUTS: { N ORDER OF THE ORIGINAL MATRIX A. { NDIM DECLARED DIMENSION OF ARRAY A. { A CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS { OF THE MATRIX ARE STORED IN THE COLUMNS OF A AND { THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS { ML+1 THROUGH 2*ML+MU+1 OF A. { ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). { MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). { OUTPUTS: { A AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND { THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. { IP INDEX VECTOR OF PIVOT INDICES. { IP[N] (-1)^(NUMBER OF INTERCHANGES) OR O (^ for Power). { IER = 0 if MATRIX A IS NONSINGULAR, OR = K if FOUND TO BE { SINGULAR AT STAGE K. { USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM. { DETERM[A] = IP[N]*A[MD,1]*A[MD,2]*...*A[MD,N] WITH MD=ML+MU+1. { if IP[N]=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO. { { REFERENCE.. { THIS IS A MODifICATION OF { C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, { C.A.C.M. 15 [1972], P. 274. {----------------------------------------------------------------------*/ *IER = 0; IP[N] = 1; MD = ML + MU + 1; MD1 = MD + 1; JU = 0; if (ML == 0) goto e70; if (N == 1) goto e70; if (N < MU+2) goto e7; for (J=MU+2; J<=N; J++) for (I=1; I<=ML; I++) A[I][J] = ZERO; e7: NM1 = N - 1; for (K=1; K<=NM1; K++) { KP1 = K + 1; M = MD; MDL = IMIN(ML,N-K) + MD; for (I=MD1; I<=MDL; I++) if (ABS(A[I][K]) > ABS(A[M][K])) M = I; IP[K] = M + K - MD; T = A[M][K]; if (M == MD) goto e20; IP[N] = -IP[N]; A[M][K] = A[MD][K]; A[MD][K] = T; e20: if (T == ZERO) goto e80; T = ONE / T; for (I=MD1; I<=MDL; I++) A[I][K] = -A[I][K]*T; JU = IMIN(IMAX(JU,MU+IP[K]),N); MM = MD; if (JU < KP1) goto e60; for (J = KP1; J<=JU; J++) { M = M - 1; MM = MM - 1; T = A[M][J]; if (M == MM) goto e30; A[M][J] = A[MM][J]; A[MM][J] = T; e30: if (T == ZERO) goto e45; JK = J - K; for (I=MD1; I<=MDL; I++) { IJK = I - JK; A[IJK][J] += A[I][K]*T; } e45: ;} e60: ;} e70: K = N; if (A[MD][N] == ZERO) goto e80; return; e80: *IER = K; IP[N] = 0; } //DECB() void SOLB(int N, int NDIM, MAT5 A, int ML, int MU, double *B, int *IP) { // Labels: e25, e50 double T; int IMD,K,KB,KMD,LM,MD,M,MD1,MDL,MDM,NM1; /*------------------------------------------------------------------ { SOLUTION OF LINEAR SYSTEM, A*X = B. { INPUTS: { N ORDER OF MATRIX A. { NDIM DECLARED DIMENSION OF ARRAY A. { A TRIANGULARIZED MATRIX OBTAINED FROM DECB. { ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). { MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). { B RIGHT HAND SIDE VECTOR. { IP PIVOT VECTOR OBTAINED FROM DECB. { DO NOT USE if DECB HAS SET IER <> 0. { OUTPUT: { B SOLUTION VECTOR, X. {------------------------------------------------------------------*/ MD = ML + MU + 1; MD1 = MD + 1; MDM = MD - 1; NM1 = N - 1; if (ML == 0) goto e25; if (N == 1) goto e50; for (K=1; K<=NM1; K++) { M = IP[K]; T = B[M]; B[M] = B[K]; B[K] = T; MDL = IMIN(ML,N-K) + MD; for (I=MD1; I<=MDL; I++) { IMD = I + K - MD; B[IMD] += A[I][K]*T; } } e25: for (KB=1; K<=NM1; K++) { K = N + 1 - KB; B[K] /= A[MD][K]; T = -B[K]; KMD = MD - K; LM = IMAX(1,KMD+1); for (I=LM; I<=MDM; I++) { IMD = I - KMD; B[IMD] += A[I][K]*T; } } e50: B[1] /= A[MD][1]; } void main() { // Initialize parameters (see ROS4) N=5; //DIMENSION OF THE SYSTEM (N=5 for Example #2) IFCN=1; //FCN(N,X,Y,F) MAY DEPEND ON X X=ZERO; //INITIAL X-VALUE XEND=1.5; //FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) H=0.001; //INITIAL STEP SIZE GUESS (0.01 FOR #2) RTOL=1e-8; //RELATIVE ERROR TOLERANCE (HERE SCALAR) ATOL=1e-10; //ABSOLUTE ERROR TOLERANCE (HERE SCALAR) ITOL=0; //BOTH RTOL AND ATOL ARE SCALARS IJAC=0; //JACOBIAN IS COMPUTED INTERNALLY BY FINITE //DIFFERENCES, Procedure JAC IS NEVER CALLED MLJAC=N; //JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA //IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IDFX=0; //DF/DX IS COMPUTED INTERNALLY BY FINITE //DIFFERENCES, Function DFX IS NEVER CALLED IMAS=0; //M IS SUPPOSED TO BE THE IDENTITY //MATRIX, Procedure MAS IS NEVER CALLED MLMAS=N; //MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA //IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IOUT=0; //Procedure SOLOUT IS NEVER CALLED LE1=N; //if MLJAC=N (FULL JACOBIAN) LJAC=N; //if MLJAC=N (FULL JACOBIAN) LMAS=0; //if IMAS=0 } LIWORK= N+2; //DECLARED LENGTH OF ARRAY IWORK LWORK = N*(LJAC+LMAS+LE1+8)+5; //DECLARED LENGTH OF ARRAY LWORK //dynamic allocations vmblock = vminit(); // initialize storage Y = (REAL *) vmalloc(vmblock, VEKTOR, NMX, 0); WORK = (REAL *) vmalloc(vmblock, VEKTOR, LWORK, 0); IWORK = (int *) vmalloc(vmblock, VEKTOR, LIWORK, 0); if (!vmcomplete(vmblock)) // out of memory? printf(" Main: out of memory.\n"); for (I=1; I<=NMX; I++) { WORK[I]=ZERO; //This triggers default values (see ROS4) IWORK[I]=0; } Y[1]=ONE; //INITIAL VALUES FOR Y Y[2]=ONE; //In #2, Y[1] = Y[2] = ... = Y[5] = 1.0 Y[3]=ONE; Y[4]=ONE; Y[5]=ONE; // call Rosenbrock Procedure with appropriate parameters // (here, FCN has been removed from input parameters). ROS4(N,IFCN,&X,Y,&XEND, &H, RTOL,ATOL,ITOL, IJAC,MLJAC,MUJAC,IDFX, IMAS,MLMAS,MUMAS, IOUT,WORK,LWORK,IWORK,LIWORK, &IDID); //print results printf("\n SOLUTION AT X = %f\n", X); for (I=1; I<=N; I++) printf(" Y(%d) = %f\n", I, Y[I]); printf("\n LAST STEP SIZE = %f\n", H); printf(" ERROR CODE = %d\n\n", IDID); vmfree(vmblock); //free memory } // end of file tros4.cpp