#include #ifndef MAX_DIM #define MAX_DIM 10 #endif int CalcExp // computes a matrix exponential // returns 0 only when the overflow protection... // ...stops the computation // else returns 1 (double a [MAX_DIM+1][MAX_DIM+1], // input matrix double e [MAX_DIM+1][MAX_DIM+1], // output matrix exponential int dim); // number of coordinates int SubStep // performs an integration sub-step // returns 0 only when the overflow protection... // ...stops the computation // else returns 1 (void (funct) // computes the forcing function... // ...and its Jacobian matrix (double [], // input coordinates int, // number of coordinates double [MAX_DIM][MAX_DIM], // Jacobian matrix double []), // vector forcing function int dim, // dimension of the system double h, // time interval double x_ []); // input and output solution int IntegrationStep // performs an integration step // the return value is 1 if the required... // ...tolerances are satisfied // else the return value is 0 (void (funct) // computes the forcing function // and its Jacobian matrix (double [], // input coordinates int, // number of coordinates double [MAX_DIM][MAX_DIM], // Jacobian matrix double []), // vector forcing function int dim, // dimension of the system double t0, // initial time double h, // time interval double h_min, // minimum value for h double h_max, // maximum value for h double tol_rel, // tolerance for the relative error double tol_abs, // tolerance for the absolute error double x_ [], // input/output solution vector FILE* res); // pointer to output file int TestTolerance // tests the numerical equality (within... // ... tolerances) of two vectors // returns 1 if at least one tolerance is satisfied // returns 0 if both tolerances are not satisfied (double xi_ [], // first vector double x_ [], // second vector int dim, // dimension of the vectors double tol_rel, // tolerance for the relative error double tol_abs); // tolerance for the absolute error void PrintResult (int dim, // dimension of the solution vector double t, // current time double x_ [], // current solution vector FILE* res); // pointer to output file int CalcExp (double a [MAX_DIM+1][MAX_DIM+1], double e [MAX_DIM+1][MAX_DIM+1], int dim) // input, output and return value are explained above { double p [MAX_DIM+2][MAX_DIM+2]; double q [MAX_DIM+2][MAX_DIM+2]; double w [MAX_DIM+2][MAX_DIM+2]; double d [MAX_DIM+2]; double r [MAX_DIM+1]; double row_sum_vec [MAX_DIM+1], col_sum_vec [MAX_DIM+1]; double bin_fac, sum, abs_sum, col_sum, row_sum, resc_fac, matrix_norm; int b; int i, j, k, l, m; int m_max, r_or_c; // tunable parameters: double ROUND_FAC = 4.0; // rounding factor double PREC = 1.0e-15; // machine precision double MIN_RESC_FAC = 2.0; // stopping criterium for rescaling iterations double UP_BOUND = 1.0e10; // overflow protection threshold double LO_BOUND = 1.0e-10; // underflow protection threshold int MAX_RESC_NUM = 100; // maximum number of one-dimensional... // ...minimisations for rescaling for (m = 0; m < dim; m++) // initialises the scale factors r [m] = 1.0; for (i = 0; i < MAX_RESC_NUM; i++) // iterations to update... // ... the rescaling matrix { for (m = 0; m < dim; m++) // computes row and column norms { col_sum = LO_BOUND; row_sum = LO_BOUND; for (k = 0; k < dim; k++) { if (k != m) { row_sum += a [m][k] * a [m][k]; col_sum += a [k][m] * a [k][m]; } } row_sum_vec [m] = row_sum; col_sum_vec [m] = col_sum; } for (resc_fac = 0.0, m = 0; m < dim; m++) // find which element of... // ...the rescaling matrix... // ... is to be updated { if (resc_fac < row_sum_vec [m] / col_sum_vec [m]) { resc_fac = row_sum_vec [m] / col_sum_vec [m]; m_max = m; r_or_c = 0; } if (resc_fac < col_sum_vec [m] / row_sum_vec [m]) { resc_fac = col_sum_vec [m] / row_sum_vec [m]; m_max = m; r_or_c = 1; } } resc_fac = sqrt (sqrt (resc_fac)); // updates... if (!r_or_c) resc_fac = 1.0 / resc_fac; r [m_max] *= resc_fac; // ...rescaling for (k = 0; k < dim; k++) // performs rescaling of matrix A { if (k != m_max) { a [m_max][k] *= resc_fac; a [k][m_max] /= resc_fac; } } if (resc_fac < MIN_RESC_FAC && resc_fac > 1.0 / MIN_RESC_FAC) break; // stops rescaling iterations } for (matrix_norm = 0.0, m = 0; m < dim; m++) // computes matrix norm { for (abs_sum = 0.0, k = 0; k < dim; k++) abs_sum += fabs (a [m][k]); if (abs_sum > matrix_norm) matrix_norm = abs_sum; } matrix_norm = matrix_norm * matrix_norm * 4.0 / PREC; for (b = 0, bin_fac = 1.0; matrix_norm > 1.0 || b <= 2; b++) // computes b... // ...and 2^(-b) { bin_fac /= 2.0; matrix_norm /= 2.0; } for (m = 0; m < dim; m++) // initialises matrix P for (k = 0; k < dim; k++) p [m][k] = a [m][k] * bin_fac; for (m = 0; m < dim; m++) // computes the last row of matrix P { for (sum = 0.0, k = 0; k < dim; k++) sum += p [k][m]; p [dim][m] = -sum; } for (m = 0; m < dim; m++) // sets the last column of matrix P p [m][dim] = 0.0; for (m = 0; m < dim; m++) // adds 1 to the main diagonal of matrix P p [m][m] += 1.0; p [dim][dim] = 1.0; for (j = 0; j < b; j++) // loop on squaring { for (i = 0; i <= dim; i++) // computes... for (k = 0; k <= dim; k++) { q [i][k] = 0.0; for (l = 0; l <= dim; l++) q [i][k] += (p [i][l] * p [l][k]); } // ...a squaring for (m = 0; m <= dim; m++) // computes weights { for (sum = 0.0, k = 0; k <= dim; k++) sum += q [k][m] * q [k][m]; if (sum > UP_BOUND) return 0; // overflow protection for (l = 0; l <= dim; l++) w [l][m] = (q [l][m] * q [l][m]) / sum; } for (m = 0; m <= dim; m++) // computes column errors { for (sum = 0.0, abs_sum = 0.0, k = 0; k <= dim; k++) { sum += q [k][m]; abs_sum += fabs (q [k][m]); } abs_sum *= ROUND_FAC; d [m] = 1.0 - sum; d [m] += abs_sum; d [m] -= abs_sum; } for (m = 0; m <= dim; m++) // smearing of column errors for (l = 0; l <= dim; l++) p [l][m] = q [l][m] + w [l][m] * d [m]; } // end of squaring loop for (m = 0; m < dim; m++) // sets output matrix E for (l = 0; l < dim; l++) e [l][m] = p [l][m]; for (m = 0; m < dim; m++) // de-rescaling of risult... for (k = 0; k < dim; k++) { if (k != m) { e [m][k] /= r [m]; e [k][m] *= r [m]; } } // ...calculating E = R^(-1) * E * R return 1; } int SubStep ( void (funct) (double [], int, double [MAX_DIM][MAX_DIM], double []), int dim, double h, double x_ [] ) // input, output and return value are explained above { double M [MAX_DIM+1][MAX_DIM+1]; double R [MAX_DIM+1][MAX_DIM+1]; double J [MAX_DIM][MAX_DIM]; double F [MAX_DIM]; int i, j; (funct) (x_, dim, J, F); // computes the forcing function... // ... and its Jacobian matrix for (i = 0; i < dim; i++) // computes... { for (j = 0; j < dim; j++) M [i][j] = h * J [i][j]; M [i][dim] = h * F [i]; M [dim][i] = 0; } M [dim][dim] = 0; // ...the augmented matrix hA // compute the matrix exponential if (!CalcExp (M, R, dim + 1)) return 0; for (i = 0; i < dim; i++) // extract the solution vectors x_ [i] += R [i][dim]; return 1; } int TestTolerance (double xi_ [], double x_ [], int dim, double tol_rel, double tol_abs) // input, output and return value are explained above { int i; double diff, med; for (i = 0; i < dim; i++) { diff = fabs (xi_ [i] - x_ [i]); med = fabs (xi_ [i] + x_ [i]) / 2; if (diff > tol_abs && diff > tol_rel * med) return 0; } return 1; } int IntegrationStep (void (funct) (double [], int, double [MAX_DIM][MAX_DIM], double []), int dim, double t0, double h, double h_min, double h_max, double tol_rel, double tol_abs, double x_ [], FILE* res) // input, output and return value are explained above { double xi_ [MAX_DIM]; double xm_ [MAX_DIM]; double xf_ [MAX_DIM]; int i; int I1, I2; if (dim > MAX_DIM) // check dimension { printf ("ERROR: dimension incorrect. dim > MAX_DIM\n"); abort (); } if (h > h_max) // if step too long:... // ...split it into two consecutive half-length steps { I1 = IntegrationStep (funct, dim, t0, h/2, h_min, h_max, tol_rel, tol_abs, x_, res); I2 = IntegrationStep (funct, dim, t0+h/2, h/2, h_min, h_max, tol_rel, tol_abs, x_, res); return (I1 && I2); } for (i = 0; i < dim; i++) xi_ [i] = x_ [i]; if (SubStep (funct, dim, 2 * h, x_)) // attempts a double-length step { for (i = 0; i < dim; i++) { xf_ [i] = x_ [i]; x_ [i] = xi_ [i]; } if (SubStep (funct, dim, h, x_)) // attempts a single-length step { for (i = 0; i < dim; i++) xm_ [i] = x_ [i]; if (SubStep (funct, dim, h, x_)) // attempts a 2nd single-length step { // if tolerance test OK, accept first single-length step... if (TestTolerance (xf_, x_, dim, tol_rel, tol_abs)) { for (i = 0; i < dim; i++) x_ [i] = xm_ [i]; PrintResult (dim, t0+h, x_, res); return 1; } } } } // ... else attempt, if possible, two consecutive half-length steps if (h < h_min) return 0; for (i = 0; i < dim; i++) x_ [i] = xi_ [i]; I1 = IntegrationStep (funct, dim, t0, h/2, h_min, h_max, tol_rel, tol_abs, x_, res); I2 = IntegrationStep (funct, dim, t0+h/2, h/2, h_min, h_max, tol_rel, tol_abs, x_, res); return (I1 && I2); }