mirror of
https://github.com/libretro/cpp-cheat.git
synced 2025-04-07 05:21:55 +00:00
96 lines
2.0 KiB
C
96 lines
2.0 KiB
C
#include "common.h"
|
|
|
|
/*
|
|
equation:
|
|
|
|
$Y'(t) = f(Y(t),t)$
|
|
|
|
where
|
|
|
|
$Y(t) = (y(t), y'(t), y''(t), ...)$
|
|
*/
|
|
int ode_func(double t, const double y[], double f[], void *params) {
|
|
double mu = *((double*)params);
|
|
f[0] = y[1];
|
|
f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
|
|
return GSL_SUCCESS;
|
|
}
|
|
|
|
/*
|
|
Some methods ask for you to give the jacobian of it too
|
|
for better performance.
|
|
|
|
You have to calculate this by hand.
|
|
*/
|
|
int ode_jac(double t, const double y[], double *dfdy, double dfdt[], void *params) {
|
|
double mu = *((double*)params);
|
|
gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
|
|
gsl_matrix* m = &dfdy_mat.matrix;
|
|
gsl_matrix_set(m, 0, 0, 0.0);
|
|
gsl_matrix_set(m, 0, 1, 1.0);
|
|
gsl_matrix_set(m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
|
|
gsl_matrix_set(m, 1, 1, -mu*(y[0]*y[0] - 1.0));
|
|
dfdt[0] = 0.0;
|
|
dfdt[1] = 0.0;
|
|
return GSL_SUCCESS;
|
|
}
|
|
|
|
int main(void) {
|
|
double mu = 1.0;
|
|
gsl_odeiv2_system sys = { ode_func, ode_jac, 2, &mu };
|
|
gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(
|
|
&sys,
|
|
gsl_odeiv2_step_rk8pd,
|
|
1e-6,
|
|
1e-6,
|
|
0.0
|
|
);
|
|
int i;
|
|
double t = 0.0;
|
|
double t1 = 100.0;
|
|
double y[2] = {1.0, 0.0};
|
|
const int nPts = 1000;
|
|
double dt = t1 / nPts;
|
|
FILE* fp = fopen(__FILE__ "out.tmp", "w");
|
|
double datax[nPts];
|
|
double datay[nPts];
|
|
|
|
for (i = 0; i < nPts; i++) {
|
|
double ti = i * dt;
|
|
int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
|
|
if (status != GSL_SUCCESS) {
|
|
fprintf(stderr, "error, return value=%d\n", status);
|
|
break;
|
|
}
|
|
|
|
/* Get output. */
|
|
fprintf(fp, "%.5e %.5e %.5e\n", t, y[0], y[1]);
|
|
printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
|
|
datax[i] = y[0];
|
|
datay[i] = y[1];
|
|
}
|
|
|
|
/* Cleanup. */
|
|
fclose(fp);
|
|
gsl_odeiv2_driver_free(d);
|
|
|
|
/* Plot. */
|
|
plsdev("xwin");
|
|
plinit();
|
|
plenv(
|
|
gsl_stats_min(datax, 1, nPts),
|
|
gsl_stats_max(datax, 1, nPts),
|
|
gsl_stats_min(datay, 1, nPts),
|
|
gsl_stats_max(datay, 1, nPts),
|
|
0,
|
|
1
|
|
);
|
|
plstring(nPts, datax, datay, "*");
|
|
plend();
|
|
|
|
/* #gnuplot */
|
|
/*system("gnuplot -e 'p \"ode.tmp\" u 2:3 w lines; pause -1'");*/
|
|
|
|
return EXIT_SUCCESS;
|
|
}
|