//************************************************************************ // Harmonical oscillator with the Heun's method // x'=y // y'=-x, x=x(t), y=y(t), T>t>0, x(0)=0, y(0)=1 //************************************************************************ /* compile with gcc -o test_heun test_heun.c -L/usr/local/dislin -ldislnc -lm -lX11 */ #include #include #include "/usr/local/dislin/dislin.h" #define NX 25 #define NY 25 // int step = 0; double max_steps = 40.0*3.1415926; // final time // double x = 0.0; //initial conditions double y = 1.0; // double h = 0.05; // time step double L=6.0; // double m_0x = 0.0; double m_1x = 0.0; double m_0y = 0.0; double m_1y = 0.0; // float xv[NX][NY], yv[NX][NY], xp[NX], yp[NY]; //define functions for r.h.s. double funct_x(double x, double y) { return y; } double funct_y(double x, double y) { return -x; } main () { int i, j, iret , iclr; double xstep, ystep; // define grid for a vector plot xstep = 2*L / (NX - 1);; ystep = 2*L / (NY - 1);; for (i = 0; i < NX; i++) { xp[i] = (float) (-L + i * xstep); for (j = 0; j < NY; j++) { yp[j] = (float) (-L + j *ystep); xv[i][j] = yp[j]; yv[i][j] = -xp[i]; } } winsiz (600, 600); page (2600, 2600); sclmod ("full"); scrmod ("revers"); metafl ("xwin"); disini (); pagera (); hwfont (); name ("x", "x"); name ("y", "y"); axspos (350, 2300); axslen (2000, 2000); titlin ("Phase diagram", 4); graf (-L, L, -L, 1., -L, L, -L, 1.); height (50); title (); setrgb(0,0,1); // set blue color vecmat ((float *) xv, (float *) yv, NX, NY, xp, yp, -1); // plot vector field //----------------------------------------------------------------- // Heun's method for ( step = 0; step < max_steps; step++ ) { m_0x = h * funct_x(x, y); m_0y = h * funct_y(x, y); m_1x = h * funct_x(x + m_0x, y + m_0y); m_1y = h * funct_y(x + m_0x, y + m_0y); x += (m_0x + m_1x) / 2.0; y += (m_0y + m_1y) / 2.0; setrgb(0,0,0); hsymbl(10); rlsymb(21,x,y); } //----------------------------------------------------------------- /*stmmod ("rk4", "integration"); stmmod ("on", "close"); color ("blue"); iclr = intrgb (0.0, 0.0, 0.0); vecclr (iclr); stream ((float *) xv, (float *) yv, NX, NY, xp, yp, NULL, NULL, 0); //plot the stream */ disfin (); }