/*
Programmazione II - Esempio di codice
<![CDATA[
Example code: Runge-Kutta integration of a 1-dim dynamic system
with adaptive stepsize calculation.
(20000816 prelz@mi.infn.it)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "graph_helper.h"
#ifndef TRUE
#define TRUE 1
#endif
#ifndef FALSE
#define FALSE 0
#endif
struct vel_pass
{
double t0;
double v0;
} ;
double fourth_order_rk_step(double t0, double dt, double q0, double dq0,
void *user_data,
double (*deriv)(double t, double q, void *user_data) );
double velocity_deriv(double t, double q, void *us_pointer);
double oscillator_force (double t, double v, void *us_pointer);
double duffing_force (double t, double v, void *us_pointer);
#define MAX_REPEATS 5
double (*force)(double t, double q, void *user_data);
int
main (int argc, char *argv[])
{
graph *mygr;
double xcur, vcur, vnew, xnew, v1, v2, tcur, tstep;
double req_trunc_error=5E-7;
double est_trunc_error;
char dumb[8];
int i;
struct vel_pass velocity;
mygr = graph_create_new(NULL,"Phase diagram",500,500,"x","v");
if (mygr == NULL)
{
fprintf(stderr,"Cannot create graph.\n");
exit(1);
}
graph_set_scale(mygr,-2.5,2.5,-1.5,1.5);
graph_clear(mygr);
force = oscillator_force;
for (i=0; i<MAX_REPEATS; i++)
{
graph_pick_values(mygr, &xcur, &vcur);
graph_clear_values(mygr);
/* Initial step guess */
tstep = 0.1;
for (tcur = 0; tcur < 50; )
{
velocity.t0 = tcur;
velocity.v0 = vcur;
/* Test whether this stepsize is any good */
vnew = fourth_order_rk_step(tcur, tstep, vcur,
(*force)(tcur, vcur, &xcur),
&xcur,
(*force));
v1 = vnew;
v2 = fourth_order_rk_step(tcur, tstep*2, vcur,
(*force)(tcur, vcur, &xcur),
&xcur,
(*force));
xnew = fourth_order_rk_step(tcur, tstep, xcur, vcur,
(void *)(&velocity),
velocity_deriv);
v1 = fourth_order_rk_step(tcur, tstep, v1,
(*force)(tcur, v1, &xnew),
&xnew,
(*force));
est_trunc_error = fabs(v1 - v2);
if (est_trunc_error > req_trunc_error)
{
/* Reduce step and waste this cycle */
tstep *= pow(req_trunc_error/est_trunc_error,0.2);
continue;
}
/* Adapt step and continue */
tcur += tstep;
xcur = xnew;
vcur = vnew;
tstep *= pow(req_trunc_error/est_trunc_error,0.2);
graph_add_and_plot(mygr, xcur, vcur, "green4");
}
printf("Final step: %g\n",tstep);
printf("Final trunc error: %g\n",est_trunc_error);
graph_plot(mygr, "red");
}
printf("Press <Enter> to finish\n");
fgets(dumb, sizeof(dumb), stdin);
graph_free(mygr);
XHLP_close_display();
exit(0);
}
double
fourth_order_rk_step(double t0, double dt, double q0, double dq0,
void *user_data,
double (*deriv)(double t, double q, void *user_data) )
{
double retval;
double k1,k2,k3,k4;
retval = q0;
k1 = dq0 * dt;
retval += k1/6;
k2 = deriv(t0 + dt/2, q0 + k1/2, user_data) * dt;
retval += k2/3;
k3 = deriv(t0 + dt/2, q0 + k2/2, user_data) * dt;
retval += k3/3;
k4 = deriv(t0 + dt, q0 + k3, user_data) * dt;
retval += k4/6;
if (retval == q0)
{
fprintf(stderr,"Error in fourth_order_rk_step: step is too small!\n");
}
return(retval);
}
double
velocity_deriv(double t, double q, void *user_pointer)
{
double retval;
static double last_t = 0;
static double last_v = 0;
static last_t_valid = FALSE;
struct vel_pass *vel0;
vel0 = (struct vel_pass *)user_pointer;
if (vel0->t0 == t)
{
retval = vel0->v0;
}
else if (last_t_valid && t == last_t)
{
retval = last_v;
}
else
{
retval = fourth_order_rk_step(vel0->t0, (t - vel0->t0), vel0->v0,
(*force)(vel0->t0, vel0->v0, &q),
&q,
(*force));
}
last_t = t;
last_v = retval;
last_t_valid = TRUE;
return(retval);
}
double
oscillator_force (double t, double v, void *us_pointer)
{
double x;
x = *(double *)us_pointer;
return (-0.5*x);
}
double
duffing_force (double t, double v, void *us_pointer)
{
double x;
x = *(double *)us_pointer;
return (x*(1-x*x) + 0.4*cos(t) - 0.25*v);
}
]]></plaintext></span>