#include #include #define pi 3.14159265 /* Time, position, and velocity define system state */ struct state { float t, x, v; }; /* Define equation of motion - compute acceleration */ float EoM(struct state ori) { float accel; accel = -1*ori.x; return accel; } /* Numerical solutions to ODEs */ struct state update_Euler(float step, struct state ori) { /* Initialize output structure */ struct state ter; /* Increment time, position, velocity */ ter.t = ori.t + step; ter.x = ori.x + step*ori.v; ter.v = ori.v + step*EoM(ori); return ter; } struct state update_RK4(float step, struct state ori) { /* Initialize output structure */ struct state ter, k1, k2, k3, k4; /* Update state in four iterative steps */ k1.t = step; k1.x = step*ori.v; k1.v = step*EoM(ori); k2.t = step/2; k2.x = step*(ori.v + k1.v/2); k2.v = step*(EoM(ori) + EoM(k1)/2); k3.t = step/2; k3.x = step*(ori.v + k2.v/2); k3.v = step*(EoM(ori) + EoM(k2)/2); k4.t = step; k4.x = step*(ori.v + k3.v); k4.v = step*(EoM(ori) + EoM(k3)); ter.t = ori.t + step; ter.x = ori.x + k1.x/6 + k2.x/3 + k3.x/3 + k4.x/6; ter.v = ori.v + k1.v/6 + k2.v/3 + k3.v/3 + k4.v/6; return ter; } int main() { /* Initialize output file */ FILE *fp; fp = fopen("7pt2.txt", "w"); /* Set parameters */ float step = 0.01, duration = 100*pi, est_error, act_error; int Nsteps = duration/step; struct state ori; ori.t = 0; ori.x = 1; ori.v = 0; /* Initialize output */ fprintf(fp, "time\tposition\tvelocity\test_error\tact_error\n"); fprintf(fp, "%f\t%f\t%f\t%f\t%f\n", ori.t, ori.x, ori.v, est_error, act_error); /* Compute trajectory */ struct state ter_full, ter_half1, ter_half2; int i; for(i = 1; i < Nsteps; ++i) { switch (2) { case 1: /* Euler method */ /* Advance full step */ ter_full = update_Euler(step, ori); /* To estimate the local error, advance two half steps */ ter_half1 = update_Euler(step/2, ori); ter_half2 = update_Euler(step/2, ter_half1); /* Estimate local error */ est_error = (ter_full.x - ter_half2.x)/pow(step, 2); break; case 2: /* fixed-step fourth-order Runge-Kutta */ /* Advance full step */ ter_full = update_RK4(step, ori); /* To estimate the local error, advance two half steps */ ter_half1 = update_RK4(step/2, ori); ter_half2 = update_RK4(step/2, ter_half1); /* Estimate local error */ est_error = (ter_full.x - ter_half2.x)/pow(step, 5); break; } /* Calculate local error based on analytical solution */ act_error = ter_full.x - cos(ter_full.t); /* Output state */ fprintf(fp, "%f\t%f\t%f\t%f\t%f\n", ter_full.t, ter_full.x, ter_full.v, est_error, act_error); /* Update ori for next loop iteration */ ori = ter_full; } fclose(fp); }