Wikimore

Summary

Description	English: Solve an ordinary differential equation and output the solution to stdout: ./ode.out Plot with gnuplot: ./ode.out > ode.dat ./ode.gnuplot xdg-open ode.png To also plot the phase space on a plplot output window use an argument: ./ode.out 1 An ODE of any order can be reduced to a system of first order ODEs: therefore, GSL only offers an API to solve systems of first order ODEs. In this particular example we, solve the second order autonomous ODE equation: g : R -> R g(t) = k * g'(t) * (1 - g(t)^2) - g(t) = f(g(t), g'(t), t) which is equivalent to solving the first order system with two functions y_0 and y_1: y_0'(t) = y_1(t) y_1'(t) = y_1(t) * (1 - y_0(t)^2) - y_0(t) with: y_0(t) = g(t) y_1(t) = g'(t) This function was chosen because it quickly converges to a limit cycle, which is fun to watch. Generated with: https://github.com/cirosantilli/cpp-cheat/blob/890c065b15a350b61b0780f5e47721b05df1f2f6/gsl/ode.md
Date	21 August 2019
Source	Own work
Author	Cirosantilli2

Licensing

I, the copyright holder of this work, hereby publish it under the following license:

This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.

You are free:

to share – to copy, distribute and transmit the work
to remix – to adapt the work

Under the following conditions:

attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

Category:CC-BY-SA-4.0#Van%20der%20Pol%20oscillator%20solution%20and%20phase%20space%20generated%20with%20GNU%20Scientific%20Library.png

Category:Self-published work

C src code

ode.c using GNU GSL and PLplot libraries

/*

gcc o.c -Wall -Wextra -lm -lplplot -lgsl -lgslcblas


*/

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_statistics.h>
#include <plplot/plplot.h>






int ode_func(double t, const double y[], double f[], void *params) {
    double k = *((double*)params);
    (void)t;
    f[0] = y[1];
    f[1] = k * y[1] * (1.0 - y[0] * y[0]) - y[0];
    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 k = *((double*)params);
    (void)t;
    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, -(1.0 + 2.0 * k * y[0] * y[1]));
    gsl_matrix_set(m, 1, 1, k * (1.0 * y[0] * y[0]));
    /* Autonomous. */
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;
    return GSL_SUCCESS;
}

int main(int argc, char **argv) {
    enum Constexpr {n_points = 1000};
    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;
    /* Initial condition: f = 0 with speed 0. */
    double y[2] = {1.0, 0.0};
    double dt = t1 / n_points;
    double datax[n_points];
    double datay[n_points];

    for (i = 0; i < n_points; 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. */
        printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
        datax[i] = y[0];
        datay[i] = y[1];
    }

    /* Cleanup. */
    gsl_odeiv2_driver_free(d);

    /* Plot. */
    if (argc > 1 && argv[1][0] == '1') {
        plsdev("xwin");
        plinit();
        plenv(
            gsl_stats_min(datax, 1, n_points),
            gsl_stats_max(datax, 1, n_points),
            gsl_stats_min(datay, 1, n_points),
            gsl_stats_max(datay, 1, n_points),
            0,
            1
        );
        plstring(n_points, datax, datay, "*");
        plend();
    }

    return EXIT_SUCCESS;
}

Gnuplot code: ode.gnuplot

#!/usr/bin/env gnuplot
set terminal png size 1024, 1024
set output "ode.png"
set multiplot layout 3,1 title "ODE solution"

set title "Solution: t x g(t) "
set xlabel "t"
set ylabel "g(t)"
plot "ode.dat" using 1:2 notitle

set title "Derivative plot: t x g'(t)"
set xlabel "t"
set ylabel "g'(t)"
plot "ode.dat" using 1:3 notitle

set title "Phase space: g(t) x g(t)"
set xlabel "g'(t)"
set ylabel "g(t)"
plot "ode.dat" using 2:3 notitle

./ode.out > ode.dat
./ode.gnuplot
xdg-open ode.png

Category:Van der Pol oscillator Category:GNU Scientific Library Category:Images with C source code Category:Images with Gnuplot source code Category:PLplot

Wikimore

File:Van der Pol oscillator solution and phase space generated with GNU Scientific Library.png

Summary

Licensing

C src code