Skip to content

Instantly share code, notes, and snippets.

@Hepic Hepic/RK
Created Dec 4, 2017

Embed
What would you like to do?
// Compile: gcc ergasia1_RKF4_2nd_order.c -lm -o run
#include <stdio.h>
#include <math.h> // add -lm flag at compilation
#define max(a, b) a > b ? a : b
#define min(a, b) a < b ? a : b
double func(double t, double y) {
return 2 * (y - 1) * (y - 1); // Set the function of your preference
}
int main() {
double t0, tF, y0, y0_deriv;
double Rmax = 0.0001, ScaleMin = 0.1, ScaleMax = 4.0; // default values.
printf("Enter t0, tF, y(0), y'(0): ");
scanf("%lf%lf%lf", &t0, &tF, &y0, &y0_deriv);
double h = pow(Rmax, 1.0 / 4.0);
double Hmin = h * 0.0001;
double t = t0, y = y0;
do {
if (t + h > tF) {
h = tF - t;
}
double m1 = h * func(t, y);
double m2 = h * func(t + 1.0 / 4.0 * h, y + 1.0 / 4.0 * m1);
double m3 = h * func(t + 3.0 / 8.0 * h, y + 3.0 / 32.0 * m1 + 9.0 / 32.0 * m2);
double m4 = h * func(t + 12.0 / 13.0 * h, y + 1932.0 / 2197.0 * m1 - 7200.0 / 2197.0 * m2 + 7296.0 / 2197.0 * m3);
double m5 = h * func(t + h, y + 439.0 / 216.0 * m1 - 8.0 * m2 + 3680.0 / 513.0 * m3 - 845.0 / 4104.0 * m4);
double m6 = h * func(t + 1.0 / 2.0 * h, y - 8.0 / 27.0 * m1 + 2.0 * m2 - 3544.0 / 2565.0 * m3 + 1859.0 / 4104.0 * m4 - 11.0 / 40.0 * m5);
double error_estimate = 1.0 / 360.0 * m1 - 128.0 / 4275.0 * m2 - 2197.0 / 75240.0 * m4 + 1.0 / 50.0 * m5 + 2.0 / 55.0 * m6;
double ratio = fabs(error_estimate) / h;
if (ratio < Rmax) {
t = t + h;
y += 25.0 / 216.0 * m1 + 1408.0 / 2565.0 * m3 + 2197.0 / 4104.0 * m4 - 1.0 / 5.0 * m5;
printf("%.10lf %.10lf %.10lf\n", t, h, y);
}
double ScaleFactor = 0.84 * pow(Rmax / ratio, 1.0 / 4.0);
ScaleFactor = max(ScaleFactor, ScaleMin);
ScaleFactor = min(ScaleFactor, ScaleMax);
h *= ScaleFactor;
} while (t != tF && h >= Hmin);
if (t == tF) {
printf("%.10lf %.10lf\n", t, y);
} else {
printf("συµβαίνει να είναι h < Hmin, κάποια ιδιαιτερότητα παρατηρείται ”κοντά” στο τρέχον t\n");
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.