12th order RK solver implemented

Pending testing
utat-ss · Jul 14, 2021 · 02d37fe · 02d37fe
1 parent afed599
commit 02d37fe
Show file tree

Hide file tree

Showing 6 changed files with 172 additions and 4 deletions.
diff --git 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
@@ -1,12 +1,15 @@
 #include <stdio.h>
 #include <stdlib.h>
 #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;
 }
diff --git 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
@@ -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
diff --git 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;
 }
diff --git 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
-Original file line number
+Diff line change
@@ -1,3 +1,5 @@
     .vscode
     *.pyc
     *.exe
+    *.csv
+    *.so