From 02d37fe860689bbbe449426e9d4c80d9408ffab4 Mon Sep 17 00:00:00 2001 From: Mingde Yin Date: Wed, 14 Jul 2021 18:53:37 -0400 Subject: [PATCH] 12th order RK solver implemented Pending testing --- .gitignore | 2 + c_code/main.c | 7 ++- c_code/orbit.c | 34 ++++++++++++- c_code/orbit.h | 2 + c_code/solvers.c | 130 +++++++++++++++++++++++++++++++++++++++++++++++ c_code/solvers.h | 1 + 6 files changed, 172 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index bdbfa15..738c089 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ .vscode *.pyc *.exe +*.csv +*.so diff --git a/c_code/main.c b/c_code/main.c index 1bfdb91..c1c8c8f 100644 --- a/c_code/main.c +++ b/c_code/main.c @@ -1,12 +1,15 @@ #include #include #include "orbit.h" +#include "solvers.h" int main(void) { double y0[6] = {(6378 + 500)*1e3, 0, 0, 0, 7.5e3, 0}; - double tSpan[2] = {0, 100000}; + double tSpan[2] = {0, 1000}; double tol = 1e-12; - orbit(tSpan, y0, "RK810.csv", tol); + orbitgeneral(RK1012vec, tSpan, y0, "RK1012.csv", tol); + orbitgeneral(RK810vec, tSpan, y0, "RK810.csv", tol); + return 0; } \ No newline at end of file diff --git a/c_code/orbit.c b/c_code/orbit.c index 2c1d2cf..9717d0f 100644 --- a/c_code/orbit.c +++ b/c_code/orbit.c @@ -15,11 +15,41 @@ void orbit(double tSpan[], double y0[], char output[], double ATOL) { clock_t start, end; start = clock(); - solution* result = RK810vec(combined_perturbations, tSpan, y0, ATOL); + solution* result = RK1012vec(combined_perturbations, tSpan, y0, ATOL); end = clock(); double cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC; - printf("Integration using RK810 solver completed in %.8f seconds. Steps taken: %d\n", cpu_time_used, result->n); + printf("Integration using RK1012 solver completed in %.8f seconds. Steps taken: %d\n", cpu_time_used, result->n); + + FILE *fptr = fopen(output, "w"); + + for (int i = 0; i < result->n; i++) { + fprintf(fptr, "%.15f, %.15f, %.15f, %.15f, %.15f, %.15f, %.15f\n", result->t[i], result->y[i][0], result->y[i][1], result->y[i][2], result->y[i][3], result->y[i][4], result->y[i][5]); + free(result->y[i]); + } + fclose(fptr); + + free(result->t); + free(result); +} + +void orbitgeneral(solution* (*solver)(void (*)(double, double[], double*), double[], double[], double), double tSpan[], double y0[], char output[], double ATOL) { + /* + Propagates orbit using one of the solvers. + solver: RK810 or RK1012 function pointer + tSpan: t0, tf of times. + y0: initial x y z vx vy vz state vector + output: filename of output csv data file + ATOL: tolerance + */ + clock_t start, end; + + start = clock(); + solution* result = solver(combined_perturbations, tSpan, y0, ATOL); + end = clock(); + double cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC; + + printf("Integration using solver completed in %.8f seconds. Steps taken: %d\n", cpu_time_used, result->n); FILE *fptr = fopen(output, "w"); diff --git a/c_code/orbit.h b/c_code/orbit.h index 03a8690..376d39c 100644 --- a/c_code/orbit.h +++ b/c_code/orbit.h @@ -1,6 +1,8 @@ #ifndef ORBIT_H #define ORBIT_H +#include "solvers.h" void orbit(double tSpan[], double y0[], char output[], double ATOL); +void orbitgeneral(solution* (*solver)(void (*)(double, double[], double*), double[], double[], double), double tSpan[], double y0[], char output[], double ATOL); #endif \ No newline at end of file diff --git a/c_code/solvers.c b/c_code/solvers.c index 9f5ffe2..dc50beb 100644 --- a/c_code/solvers.c +++ b/c_code/solvers.c @@ -125,5 +125,135 @@ solution* RK810vec(void (*f)(double, double[], double*), double tSpan[], double // record # of steps after finish result->n = step; + return result; +} + +solution* RK1012vec(void (*f)(double, double[], double*), double tSpan[], double y0[], double ATOL) { + /* + Integrates along time points using an RK1012 method using n timesteps on a vector of dimension m + Source: https://sce.uhcl.edu/rungekutta/rk1210.txt + + f: function for which y' = f(t, y, ARR), where ARR stores the resulting vector + tSpan: stores (t0, tf) + y0: array of length m containing initial conditions + m: dimension of vector + ord: order of method + fname: integrator filename + ATOL: tolerance + */ + + if (tSpan[1] <= tSpan[0]) return NULL; // error case + + ////// Initial estimate for n + int n = (int)((tSpan[1] - tSpan[0])/H0); // 1 step every 10 seconds + int step = 0; + + ////// Allocate result struct + solution* result = (solution*) malloc(sizeof(solution)); + result->y = malloc(n * sizeof(double*)); + result->t = malloc(n * sizeof(double)); + for (int i = 0; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double)); + + ////// Prepare Step variables + double h = H0; // initial stepsize + double h_old; // previous step variable + double t = tSpan[0]; // intial time + + // Init y0, t0 (step 0) + for (int j = 0; j < VEC_SIZE; j++) result->y[0][j] = y0[j]; + result->t[0] = t; + + // RK Step Parameters + // a coefficients (written at beta in the paper) [intermediate solution weights] + double a[25][24] ={{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-0.21604938271604937, 0.7716049382716049, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.20833333333333334, 0.0, 0.625, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.19333333333333333, 0.0, 0.22, -0.08, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.1, 0.0, 0.0, 0.4, 0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.10336447165001048, 0.0, 0.0, 0.12405309452894676, 0.4831711675610329, -0.038753024569476324, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.12403826143183333, 0.0, 0.0, 0.0, 0.21705063219795848, 0.013745579207596677, -0.06610953172676828, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.0914774894856883, 0.0, 0.0, 0.0, 0.0, -0.005443485237174697, 0.06807168016884535, 0.40839431558264105, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.08900136525025511, 0.0, 0.0, 0.0, 0.0, 0.004995282266455323, 0.397918238819829, 0.4279302107525766, -0.0865117637557827, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.06950876241349076, 0.0, 0.0, 0.0, 0.0, 0.12914694190017645, 1.530736381023113, 0.57787476112914, -0.9512947723210889, -0.40827664296563193, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.044486140329513583, 0.0, 0.0, 0.0, 0.0, -0.0038047686705696172, 0.01069550640296242, 0.020961624449990432, -0.023314602325932177, 0.0026326598106453697, 0.0031547276897702504, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.019458881511975546, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.785129491718125e-05, -4.297958590492736e-05, 1.7635898226028515e-05, 0.0653866627415027, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.2068368356642771, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.016679606710415646, -0.008795015632007103, 0.003466754553624639, -0.8612644601057177, 0.9086518820740502, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.0203926084654484, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0869469392016686, -0.019164963041014983, 0.006556291594936633, 0.09874761281274348, 0.005353646955249961, 0.3011678640109679, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0.2284104339177781, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.4987074007930252, 0.1348411683357245, -0.03874582440558342, -1.2747325747347484, 1.4391636446287717, -0.21400746796799025, 0.9582024177544303, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {2.002224776559742, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.067018099615249, 0.6239781360861395, -0.046228368550031144, -8.849732883626496, 7.7425770785085595, -0.5883585192508692, -1.1068373336238064, -0.929529037579204, 0, 0, 0, 0, 0, 0, 0, 0}, + {3.1378953341207345, 0.0, 0.0, 0.0, 0.0, 0.12914694190017645, 1.530736381023113, 0.57787476112914, 5.420882630551267, 0.2315469260348293, 0.07592929955789135, -12.372997338018651, 9.854558834647696, 0.08591114313704365, -5.652427528626439, -1.9430093524281962, -0.12835260184940453, 0, 0, 0, 0, 0, 0, 0}, + {1.3836005443219601, 0.0, 0.0, 0.0, 0.0, 0.004995282266455323, 0.397918238819829, 0.4279302107525766, -1.3029910742447577, 0.661292278669377, -0.14455977430695435, -6.965760347317982, 6.6580854323599175, -1.669973751088415, 2.064137023180353, -0.6747439626443065, -0.001156188347949395, -0.005440579086770074, 0, 0, 0, 0, 0, 0}, + {0.9512362970482877, 0.0, 0.0, 0.0, 0.21705063219795848, 0.013745579207596677, -0.06610953172676828, 0.0, 0.15228169673641445, -0.33774101835759984, -0.019282598163399577, -3.682592696968668, 3.1619787040698206, -0.3704625221068853, -0.05149742003654404, -0.0008296255321201529, 2.798010414192786e-06, 0.041860391641236026, 0.27908425509087736, 0, 0, 0, 0, 0}, + {0.10336447165001048, 0.0, 0.0, 0.12405309452894676, 0.4831711675610329, -0.038753024569476324, 0.0, -0.43831382036112243, 0.0, -0.21863663372167666, -0.031233476439471924, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.031233476439471924, 0.21863663372167666, 0.43831382036112243, 0, 0, 0, 0}, + {0.19333333333333333, 0.0, 0.22, -0.08, 0.0, 0.0, 0.0984256130499316, -0.19641088922305466, 0.0, 0.43645793049306875, 0.06526137216757211, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.06526137216757211, -0.43645793049306875, 0.19641088922305466, -0.0984256130499316, 0, 0, 0}, + {-0.21604938271604937, 0.7716049382716049, 0.0, 0.0, -0.6666666666666666, 0.0, -0.39069646929597845, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.39069646929597845, 0.6666666666666666, 0, 0}, + {0.2, 0.0, -0.1646090534979424, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1646090534979424, 0}, + {1.4717872488111041, 0.7875, 0.4212962962962963, 0.0, 0.2916666666666667, 0.0, 0.34860071762832956, 0.22949954476899484, 5.79046485790482, 0.4185875118565069, 0.307039880222474, -4.687009053506033, 3.1357166559380225, 1.4013482971096571, -5.52931101439499, -0.8531382355080633, 0.10357578037361014, -0.14047441695060095, -0.4185875118565069, -0.22949954476899484, -0.34860071762832956, -0.2916666666666667, -0.4212962962962963, -0.7875}}; + // c coefficients (written as a in the paper) {intermediate timestep values} + double c[25] = { 0.0, 0.2, 0.5555555555555556, 0.8333333333333334, 0.3333333333333333, 1.0, 0.6718357091705138, 0.2887249411106202, 0.5625, 0.8333333333333334, 0.9476954311791993, 0.054811287686380265, 0.08488805186071653, 0.2655756032646429, 0.5, 0.7344243967353571, 0.9151119481392834, 0.9476954311791993, 0.8333333333333334, 0.2887249411106202, 0.6718357091705138, 0.3333333333333333, 0.5555555555555556, 0.2, 1.0}; + + // b coefficeints (written as c in the paper) [solution weights] + double b[25] = {0.023809523809523808, 0.0234375, 0.03125, 0.0, 0.041666666666666664, 0.0, 0.05, 0.05, 0.0, 0.1, 0.07142857142857142, 0.0, 0.13841302368078298, 0.2158726906049313, 0.2438095238095238, 0.2158726906049313, 0.13841302368078298, -0.07142857142857142, -0.1, -0.05, -0.05, -0.041666666666666664, -0.03125, -0.0234375, 0.023809523809523808}; + + // Create Step Variables + double k[25][VEC_SIZE]; // k values + double y_i[VEC_SIZE]; // function inputs at each step + double y_curr[VEC_SIZE]; // current solution array, and also used to update itself for next step + + double err; // magnitude of error + + /////// INTEGRATION STEPS + while (t-h < tSpan[1]) { + for (int j = 0; j < VEC_SIZE; j++) y_curr[j] = result->y[step][j]; // get current y (for vectorization purposes) + + f(t, y_curr, k[0]); // RK Stage 0 + + // Perform all RK Stages [1, 25) + for (int r = 1; r < 25; r++) { + //// Prepare input vector + for (int j = 0; j < VEC_SIZE; j++) { + y_i[j] = y_curr[j]; // take current sol + for (int w = 0; w < r; w++) y_i[j] += h * k[w][j] * a[r][w]; // Add previous steps + } + f(t + h * c[r], y_i, k[r]); // evaluate next k + } + + // Calculate error using error estimate formula from paper + // take abs value of all errors as convervative estimate + err = 0; + for (int j = 0; j < VEC_SIZE; j++) err += fabs(49.0/640.0 * h * (k[1][j] - k[23][j])); + + //// Step size adjustment + // Determine new step size + h_old = h; + h = 0.9 * h * fmin(fmax(pow(ATOL/err, 1.0/11.0), 0.1), 10); // next step size (based on eleventh order local error) + // Min/max corrections to stop step from changing too quickly + + if (err < ATOL) { + /// step within tolerance, append solution + // check array size and increase size if needed + if ((step + 1) >= n) { + n *= 2; + result->y = realloc(result->y, n * sizeof(double*)); + result->t = realloc(result->t, n * sizeof(double)); + for (int i = step + 1; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double)); + } + + // Append to result + for (int j = 0; j < VEC_SIZE; j++) { + for (int r = 0; r < 25; r++) y_curr[j] += h_old * b[r] * k[r][j]; // Add all weights + } + for (int j = 0; j < VEC_SIZE; j++) result->y[step+1][j] = y_curr[j]; // write (separated for vectorization) + + step++; // advance step + t += h_old; // advance time + result->t[step] = t; // record time + } + // Otherwise, retry step + } + // record # of steps after finish + result->n = step; + return result; } \ No newline at end of file diff --git a/c_code/solvers.h b/c_code/solvers.h index d29b340..cd5017c 100644 --- a/c_code/solvers.h +++ b/c_code/solvers.h @@ -8,5 +8,6 @@ typedef struct AdaptiveSolution { } solution; solution* RK810vec(void (*f)(double, double[], double*), double tSpan[], double y0[], double ATOL); +solution* RK1012vec(void (*f)(double, double[], double*), double tSpan[], double y0[], double ATOL); #endif \ No newline at end of file