/***************************************************************** * bounce.c: mass-and-spring bounce simulator * * * * See http://www.jw-stumpel/bounce.html about what this * * all means. * * * * The number of masses (and springs), N, is entered as a * * parameter on the command line. Without a parameter, its * * default value is 2. * * * * Compile with * * cc -o bounce bounce.c -lgsl -lgslcblas -lm * * * * You need to have the GNU Scientific Library (gsl) installed * * (Debian package libgsl0-dev, but available on many * * other platforms as well). * * * * The output is a gnuplot file. It is human-readable to show * * the solutions, but it can also be used as input to gnuplot: * * * * ./bounce [N] >plot.txt * * gnuplot plot.txt * * * * For some reason it WILL NOT WORK if you just say * * * * ./bounce | gnuplot * * * * Please edit the output file by hand if you want to * * generate pictures in other media (png, PostScript). * * * * Some extra information will also be displayed on the * * terminal (stderr). * * * * (c) J.W. Stumpel, 14 October 2006 (jstumpel@planet.nl) * * (to whom questions should be addressed). * * * * Licence: GPL. Well, it's free to use and change. Teachers * * who use this in their classes are encouraged to send me an * * e-mail. * * * * Acknowledgements: thanks to the gsl and gnuplot people. * * Thanks to Wen-Chyuan Yueh for clarification about his * * formulae. * *****************************************************************/ #include #include #include #include #include int N = 2; gsl_matrix *A, *evec; gsl_vector *omega, *weights; double position(double t, int n) /* returns position of ball n at time t*/ { int j; double a, b, c, result = 0.0; for (j = 0; j < N; j++) { a = gsl_vector_get(weights, j); b = gsl_matrix_get(evec, n, j); c = gsl_vector_get(omega, j); result += a * b * sin(c * t); } return result; } double speed(double t, int n) /* returns speed of ball n at time t */ { int j; double a, b, c, result = 0.0; for (j = 0; j < N; j++) { a = gsl_vector_get(weights, j); b = gsl_matrix_get(evec, n, j); c = gsl_vector_get(omega, j); result += a * b * c * cos(c * t); } return result; } double findroot(void) /* finds the moment when mass 0 starts moving up again, (= the duration of the impact) by Newton-Raphson */ { int j = 0; double test = 2.0, diff; /* 2.0 is initial "guess" value */ while (1) { diff = position(test, 0) / speed(test, 0); if (fabs(diff) < 0.0000000001) /* enough precision? */ return test; test -= diff; if (++j > 20) { fprintf(stderr, "root does not converge!\n"); exit(1); } } } int main(int argc, char *argv[]) { int i, j; int s; double T, avspeed = 0, Epot = 0.0, Ekin = 0.0, Etot = 0.0, a, b; if (argc == 2) sscanf(argv[1], "%d", &N); if (!(N >= 1)) N = 2; /* Set up the initial 'A' matrix */ A = gsl_matrix_alloc(N, N); gsl_matrix_set_zero(A); gsl_matrix_set(A, N - 1, N - 1, 1.0); if (N > 1) { for (j = 1; j < N; j++) { gsl_matrix_set(A, j - 1, j, -1.0); gsl_matrix_set(A, j, j - 1, -1.0); } for (j = 0; j < N - 1; j++) gsl_matrix_set(A, j, j, 2.0); } /* compute eigenvalues and eigenvectors. This could be done with Yueh's formulas as well, but here we just use gsl. */ gsl_vector *eval = gsl_vector_alloc(N); evec = gsl_matrix_alloc(N, N); gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(N); gsl_eigen_symmv(A, eval, evec, w); gsl_eigen_symmv_free(w); gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC); /* determine the weights */ omega = gsl_vector_alloc(N); weights = gsl_vector_alloc(N); gsl_vector *V = gsl_vector_alloc(N); gsl_vector_set_all(V, -1.0); for (j = 0; j < N; j++) gsl_vector_set(omega, j, (double) N * sqrt(gsl_vector_get(eval, j))); gsl_matrix *U = gsl_matrix_alloc(N, N); gsl_matrix_memcpy(U, evec); /* copy; will be destroyed */ gsl_permutation *perm = gsl_permutation_alloc(N); gsl_linalg_LU_decomp(U, perm, &s); gsl_linalg_LU_solve(U, perm, V, weights); gsl_vector_div(weights, omega); gsl_permutation_free(perm); /* Determine the duration of the impact */ T = findroot(); /* calculate the energies involved */ for (j = 0; j < N; j++) { a = speed(T, j); avspeed += a; Ekin += a * a; } for (j = 0; j < N - 1; j++) { a = position(T, j); a -= position(T, j + 1); Epot += a * a; } Ekin /= N; /* mass scaling law */ Epot *= N; /* spring constant scaling law */ Etot = Ekin + Epot; avspeed /= N; /* coefficient of restitution */ fprintf(stderr, "Number of masses: %d\n", N); fprintf(stderr, "Duration of impact: %1.9f\n", T); fprintf(stderr, "Restitution coeff.: %1.9f\n", avspeed); fprintf(stderr, "Kinetic energy: %1.9f\n", Ekin); fprintf(stderr, "Potential energy: %1.9f\n", Epot); fprintf(stderr, "Total energy (must be 1): %1.9f\n", Etot); a = avspeed * avspeed; fprintf(stderr, "External energy: %1.9f\n", a); fprintf(stderr, "Internal (\"lost\") energy: %1.9f\n\n", 1.0 - a); /* The internal energy after the bounce (at moment T) is the potential energy + the kinetic energy which does not contribute to the common center-of-gravity movement */ /* Output section: prepares a gnuplot file. Some "scaling" of the output is applied to show the movements more clearly. Maybe there is an easier C interface to gnuplot, but I hardly know anything about gnuplot, so I just do it this way. */ printf("set terminal x11\n"); printf("set style line 1 lt 1 lw 1\n"); printf("set key off\n"); printf("set yrange [0:]\n"); printf("set title \"Bounce of a chain of %d %s and %d %s\"\n", N, (N > 1) ? "masses" : "mass", N, (N > 1) ? "springs" : "spring"); printf("plot [0:%1.4f] \\\n", T); for (i = 0; i < N; i++) { printf("%d+%3.1f*(\\\n", i + 1, 0.3 * N); /* here is the "scaling" */ for (j = 0; j < N; j++) printf("%+1.6f*sin(%1.6f*x)\\\n", gsl_vector_get(weights, j) * gsl_matrix_get(evec, i, j), gsl_vector_get(omega, j)); printf(") ls 1%s", (i == N - 1) ? "\n" : ",\\\n"); } printf("pause -1\n"); return 0; }