# include # include # include # include int main ( ); int i4_modp ( int i, int j ); int i4_wrap ( int ival, int ilo, int ihi ); double *lorenz96_rhs ( double t, int m, double x[] ); double r8_normal_01 ( ); double *r8vec_linspace_new ( int n, double a, double b ); double *rk4vec ( double t0, int n, double u0[], double dt, double *f ( double t, int n, double u[] ) ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: lorenz96_ode() solves the Lorenz96 ODE. Licensing: This code is distributed under the MIT license. Modified: 30 September 2025 Author: John Burkardt */ { char command_filename[] = "lorenz96_ode_commands.txt"; FILE *command_unit; char data_filename[] = "lorenz96_ode_data.txt"; FILE *data_unit; double dt; double force = 8.0; int i; int j; int n = 4; int nt = 2000; double perturb = 0.001; double *t; double t0 = 0.0; double tstop = 30.0; double *y; double *ynew; timestamp ( ); printf ( "\n" ); printf ( "lorenz96_ode():\n" ); printf ( " C version\n" ); printf ( " Compute solutions of the Lorenz96 ordinary differential equations (ODE).\n" ); printf ( " Write data to a file for use by gnuplot().\n" ); /* Initialize random number generator. */ srand48 ( time ( NULL ) ); /* Data */ dt = ( tstop - t0 ) / ( double ) ( nt ); /* Store the initial conditions in entry 0. */ t = r8vec_linspace_new ( nt + 1, t0, tstop ); y = ( double * ) malloc ( n * ( nt + 1 ) * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { y[i+0*nt] = force + perturb * r8_normal_01 ( ); } /* Compute the approximate solution at equally spaced times. */ for ( j = 0; j < nt; j++ ) { ynew = rk4vec ( t[j], n, y+j*n, dt, lorenz96_rhs ); for ( i = 0; i < n; i++ ) { y[i+(j+1)*n] = ynew[i]; } free ( ynew ); } /* Create the plot data file. */ data_unit = fopen ( data_filename, "wt" ); for ( j = 0; j <= nt; j++ ) { fprintf ( data_unit, " %14.6g %14.6g %14.6g %14.6g %14.6g\n", t[j], y[0+j*n], y[1+j*n], y[2+j*n], y[3+j*n] ); } fclose ( data_unit ); printf ( " Created data file \"%s\".\n", data_filename ); /* Create the plot command file. */ command_unit = fopen ( command_filename, "wt" ); fprintf ( command_unit, "# %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "# Usage:\n" ); fprintf ( command_unit, "# gnuplot < %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "set term png\n" ); fprintf ( command_unit, "set output 'lorenz96_time.png'\n" ); fprintf ( command_unit, "set xlabel '<--- T --->'\n" ); fprintf ( command_unit, "set ylabel '<--- Y0(T), Y1(T), Y2(T) --->'\n" ); fprintf ( command_unit, "set title 'Y0, Y1, Y2 versus Time'\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "plot '%s' using 1:2 lw 3 linecolor rgb 'blue',", data_filename ); fprintf ( command_unit, "'' using 1:3 lw 3 linecolor rgb 'red'," ); fprintf ( command_unit, "'' using 1:4 lw 3 linecolor rgb 'green'\n" ); fprintf ( command_unit, "set output 'lorenz96_3d.png'\n" ); fprintf ( command_unit, "set xlabel '<--- Y0(T) --->'\n" ); fprintf ( command_unit, "set ylabel '<--- Y1(T) --->'\n" ); fprintf ( command_unit, "set zlabel '<--- Y2(T) --->'\n" ); fprintf ( command_unit, "set title '(Y0(T),Y1(T),Y2(T)) trajectory'\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "splot '%s' using 2:3:4 lw 1 linecolor rgb 'blue'\n", data_filename ); fprintf ( command_unit, "quit\n" ); fclose ( command_unit ); printf ( " Created command file '%s'\n", command_filename ); /* Terminate. */ printf ( "\n" ); printf ( "lorenz96_ode():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ int i4_modp ( int i, int j ) /******************************************************************************/ /* Purpose: i4_modp() returns the nonnegative remainder of I4 division. Discussion: If NREM = I4_MODP ( I, J ) NMULT = ( I - NREM ) / J then I = J * NMULT + NREM where NREM is always nonnegative. The MOD function computes a result with the same sign as the quantity being divided. Thus, suppose you had an angle A, and you wanted to ensure that it was between 0 and 360. Then mod(A,360) would do, if A was positive, but if A was negative, your result would be between -360 and 0. On the other hand, I4_MODP(A,360) is between 0 and 360, always. Example: I J MOD I4_MODP I4_MODP Factorization 107 50 7 7 107 = 2 * 50 + 7 107 -50 7 7 107 = -2 * -50 + 7 -107 50 -7 43 -107 = -3 * 50 + 43 -107 -50 -7 43 -107 = 3 * -50 + 43 Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 Author: John Burkardt Input: int I, the number to be divided. int J, the number that divides I. Output: int I4_MODP, the nonnegative remainder when I is divided by J. */ { int value; if ( j == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "i4_modp(): Fatal error!\n" ); fprintf ( stderr, " i4_modp ( i, j ) called with J = %d\n", j ); exit ( 1 ); } value = i % j; if ( value < 0 ) { value = value + abs ( j ); } return value; } /******************************************************************************/ int i4_wrap ( int ival, int ilo, int ihi ) /******************************************************************************/ /* Purpose: i4_wrap() forces an I4 to lie between given limits by wrapping. Example: ILO = 4, IHI = 8 I Value -2 8 -1 4 0 5 1 6 2 7 3 8 4 4 5 5 6 6 7 7 8 8 9 4 10 5 11 6 12 7 13 8 14 4 Licensing: This code is distributed under the MIT license. Modified: 26 December 2012 Author: John Burkardt Input: int IVAL, an integer value. int ILO, IHI, the desired bounds for the integer value. Output: int I4_WRAP, a "wrapped" version of IVAL. */ { int jhi; int jlo; int value; int wide; if ( ilo < ihi ) { jlo = ilo; jhi = ihi; } else { jlo = ihi; jhi = ilo; } wide = jhi + 1 - jlo; if ( wide == 1 ) { value = jlo; } else { value = jlo + i4_modp ( ival - jlo, wide ); } return value; } /******************************************************************************/ double *lorenz96_rhs ( double t, int n, double y[] ) /******************************************************************************/ /* Purpose: lorenz96_rhs() evaluates the right hand side of the Lorenz96 ODE. Licensing: This code is distributed under the MIT license. Modified: 30 September 2025 Author: John Burkardt Input: double T, the value of the independent variable. int N, the spatial dimension. double Y[N], the values of the dependent variables at time T. Output: double DYDT[N], the values of the derivatives of the dependent variables at time T. */ { double *dydt; double force = 8.0; int i; int im1; int im2; int ip1; dydt = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { ip1 = i4_wrap ( i + 1, 0, n - 1 ); im1 = i4_wrap ( i - 1, 0, n - 1 ); im2 = i4_wrap ( i - 2, 0, n - 1 ); printf ( "im2, im1, i, ip1: %d %d %d %d\n", im2, im1, i, ip1 ); dydt[i] = ( y[ip1] - y[im2] ) * y[im1] - y[i] + force; } return dydt; } /******************************************************************************/ double r8_normal_01 ( ) /******************************************************************************/ /* Purpose: r8_normal_01() samples the standard normal probability distribution. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. The Box-Muller method is used, which is efficient, but generates two values at a time. Licensing: This code is distributed under the MIT license. Modified: 30 September 2025 Author: John Burkardt Output: double R8_NORMAL_01, a normally distributed random value. */ { double r1; double r2; const double r8_pi = 3.141592653589793; double x; r1 = drand48 ( ); r2 = drand48 ( ); x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 ); return x; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_linspace_new() creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Input: int N, the number of entries in the vector. double A, B, the first and last entries. Output: double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ double *rk4vec ( double t0, int m, double u0[], double dt, double *f ( double t, int n, double u[] ) ) /******************************************************************************/ /* Purpose: rk4vec() takes one Runge-Kutta step for a vector ODE. Discussion: It is assumed that an initial value problem, of the form du/dt = f ( t, u ) u(t0) = u0 is being solved. If the user can supply current values of t, u, a stepsize dt, and a function to evaluate the derivative, this function can compute the fourth-order Runge Kutta estimate to the solution at time t+dt. Licensing: This code is distributed under the MIT license. Modified: 09 October 2013 Author: John Burkardt Input: double T0, the current time. int M, the dimension of the space. double U0[M], the solution estimate at the current time. double DT, the time step. double *F ( double T, int M, double U[] ), a function which evaluates the derivative, or right hand side of the problem. Output: double RK4VEC[M], the fourth-order Runge-Kutta solution estimate at time T0+DT. */ { double *f0; double *f1; double *f2; double *f3; int i; double t1; double t2; double t3; double *u; double *u1; double *u2; double *u3; /* Get four sample values of the derivative. */ f0 = f ( t0, m, u0 ); t1 = t0 + dt / 2.0; u1 = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { u1[i] = u0[i] + dt * f0[i] / 2.0; } f1 = f ( t1, m, u1 ); t2 = t0 + dt / 2.0; u2 = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { u2[i] = u0[i] + dt * f1[i] / 2.0; } f2 = f ( t2, m, u2 ); t3 = t0 + dt; u3 = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { u3[i] = u0[i] + dt * f2[i]; } f3 = f ( t3, m, u3 ); /* Combine them to estimate the solution. */ u = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { u[i] = u0[i] + dt * ( f0[i] + 2.0 * f1[i] + 2.0 * f2[i] + f3[i] ) / 6.0; } /* Free memory. */ free ( f0 ); free ( f1 ); free ( f2 ); free ( f3 ); free ( u1 ); free ( u2 ); free ( u3 ); return u; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }