# include # include # include # include # include "ppc_array.h" # include "ppc_neural_pde.h" /******************************************************************************/ void Neural_Net_end ( struct Neural_Net_PDE *nn ) /******************************************************************************/ /* Purpose: Neural_Net_end() frees the memory associated with a neural network. Licensing: This code is distributed under the MIT license. Modified: 15 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. */ { free_vector ( nn->weights ); return; } /******************************************************************************/ double Neural_Net_error_vs_exact ( struct Neural_Net_PDE *nn, int m, int n ) /******************************************************************************/ /* Purpose: Neural_Net_error_vs_exact() computes the exact error. Licensing: This code is distributed under the MIT license. Modified: 16 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. Output: double value: the maximum error. */ { if ( nn->exact_sol == NULL ) { fprintf ( stderr, "Unable to calculate error; no exact solution given.\n" ); exit ( EXIT_FAILURE ); } double max_err = 0.0; double a = nn->bb_xrange[0]; double b = nn->bb_xrange[1]; double c = nn->bb_yrange[0]; double d = nn->bb_yrange[1]; for ( int i = 0; i <= m; i++ ) { double x = a + ( b - a ) / m * i; for ( int j = 0; j <= n; j++ ) { double y = c + ( d - c ) / n * j; nn->phi_func ( nn, x, y ); if ( nn->phi[0][0] < 0 ) { continue; } Neural_Net_eval ( nn, x, y ); double z = nn->N[0][0] * nn->phi[0][0]; double err = fabs ( z - nn->exact_sol ( x, y ) ); if ( max_err < err ) { max_err = err; } } } return max_err; } /******************************************************************************/ void Neural_Net_eval ( struct Neural_Net_PDE *nn, double x, double y ) /******************************************************************************/ /* Purpose: Neural_Net_eval() evaluates a neural network. Licensing: This code is distributed under the MIT license. Modified: 15 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. double x, y: the evaluation point. Output: struct Neural_Net_PDE *nn: the N[][] components have been evaluated. */ { int q = nn->q; double *u = nn->weights; double *v = nn->weights + q; double *w1 = nn->weights + 2 * q; double *w2 = nn->weights + 3 * q; for ( int i = 0; i <= 2; i++ ) { for ( int j = 0; j <= 2; j++ ) { nn->N[i][j] = 0.0; } } for ( int k = 0; k < q; k++ ) { double z = u[k] + w1[k] * x + w2[k] * y; sigmoid ( z, nn->sigma ); nn->N[0][0] = nn->N[0][0] + v[k] * nn->sigma[0]; nn->N[1][0] = nn->N[1][0] + v[k] * nn->sigma[1] * w1[k]; nn->N[0][1] = nn->N[0][1] + v[k] * nn->sigma[1] * w2[k]; nn->N[2][0] = nn->N[2][0] + v[k] * nn->sigma[2] * w1[k] * w1[k]; nn->N[1][1] = nn->N[1][1] + v[k] * nn->sigma[2] * w1[k] * w2[k]; nn->N[0][2] = nn->N[0][2] + v[k] * nn->sigma[2] * w2[k] * w2[k]; } return; } /******************************************************************************/ void Neural_Net_init ( struct Neural_Net_PDE *nn ) /******************************************************************************/ /* Purpose: Neural_Net_init() allocates and initializes a neural network. Licensing: This code is distributed under the MIT license. Modified: 15 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. Ooutput: struct Neural_Net_PDE *nn: the structure has been initialized. */ { nn->nweights = 4 * nn->q; make_vector ( nn->weights, nn->nweights ); /* Initialize the weights with random values in [-0.5,+0.5]. */ for ( int i = 0; i < nn->nweights; i++ ) { nn->weights[i] = ( double ) rand ( ) / RAND_MAX - 0.5; } return; } /******************************************************************************/ void Neural_Net_plot_with_maple ( struct Neural_Net_PDE *nn, int n, int m, char *outfile ) /******************************************************************************/ /* Purpose: Neural_Net_plot_with_maple() plots the solution over a grid. Discussion: The code writes a Maple script to plot the solution sampled on a grid of (n+1)x(m+1) points. Licensing: This code is distributed under the MIT license. Modified: 15 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. int n, m: the grid density in the x and y directions. char *outfile: the name of the script file. */ // We are in the domain if phi > -EPS. #define EPS 1e-12 { double a = nn->bb_xrange[0]; double b = nn->bb_xrange[1]; double c = nn->bb_yrange[0]; double d = nn->bb_yrange[1]; FILE *fp; fp = fopen ( outfile, "w" ); if ( fp == NULL ) { fprintf ( stderr, "unable to open file `%s' for writing\n", outfile ); return; } fprintf ( fp, "plots:-display(MESH([" ); for ( int i = 0; i <= m; i++ ) { double y = c + ( d - c ) / m * i; fprintf ( fp, "[" ); for ( int j = 0; j <= m; j++ ) { double x = a + ( b - a ) / n * j; nn->phi_func ( nn, x, y ); if ( - EPS < nn->phi[0][0] ) { Neural_Net_eval ( nn, x, y ); double z = nn->phi[0][0]*nn->N[0][0]; fprintf ( fp, "[%g,%g,%g],", x, y, z ); } else { fprintf ( fp, "[%g,%g,undefined],", x, y ); } } fprintf ( fp, "NULL]," ); } fprintf ( fp, "NULL]), labels=[x,y,z]);\n" ); fclose ( fp ); return; } /******************************************************************************/ void Neural_Net_plot_with_matlab ( struct Neural_Net_PDE *nn, int n, int m, char *outfile ) /******************************************************************************/ /* Purpose: Neural_Net_plot_with_matlab() plots the solution over a grid. Discussion: The code writes a MATLAB script to plot the solution sampled on a grid of (n+1)x(m+1) points. Licensing: This code is distributed under the MIT license. Modified: 16 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. int n, m: the grid density in the x and y directions. char *outfile: the name of the script file. */ { // names of matlab's plotting variables # define XVEC "PLOT_WITH_MATLAB_FD_XVEC" # define YVEC "PLOT_WITH_MATLAB_FD_YVEC" # define ZMAT "PLOT_WITH_MATLAB_FD_ZMAT" # define EPS 1e-12 // We are in the domain if phi > -EPS. double a = nn->bb_xrange[0]; double b = nn->bb_xrange[1]; double c = nn->bb_yrange[0]; double d = nn->bb_yrange[1]; FILE *fp; fp = fopen ( outfile, "w" ); if ( fp == NULL ) { fprintf ( stderr, "unable to open file `%s' for writing\n", outfile ); return; } fprintf ( fp, "%s = [\n", XVEC ); for ( int i = 0; i <= n; i++ ) { double x = a + ( b - a ) / n * i; fprintf ( fp, "%g ", x ); } fprintf ( fp, "\n];\n" ); fprintf ( fp, "%s = [\n", YVEC ); for ( int j = 0; j <= m; j++ ) { double y = c + ( d - c ) / m * j; fprintf ( fp, "%g ", y ); } fprintf ( fp, "\n];\n" ); fprintf ( fp, "%s = [\n", ZMAT ); for ( int i = 0; i <= m; i++ ) { double y = c + ( d - c ) / m * i; for ( int j = 0; j <= n; j++ ) { double x = a + ( b - a ) / n * j; nn->phi_func ( nn, x, y ); if ( - EPS < nn->phi[0][0] ) { Neural_Net_eval ( nn, x, y ); double z = nn->phi[0][0] * nn->N[0][0]; fprintf ( fp, "%g ", z ); } else { fprintf ( fp, "NaN "); } } fprintf ( fp, "\n" ); } fprintf ( fp, "];\n" ); fprintf ( fp, "surf(%s, %s, %s)\n", XVEC, YVEC, ZMAT ); fprintf ( fp, "xlabel('x')\n" ); fprintf ( fp, "ylabel('y')\n" ); fprintf ( fp, "zlabel('u')\n" ); fprintf ( fp, "colormap jet\n" ); fprintf ( fp, "print ( '-dpng', 'ppc_neural_pde.png' )\n" ); fclose ( fp ); return; } /******************************************************************************/ double Neural_Net_residual ( double *weights, int nweights, void *params ) /******************************************************************************/ /* Purpose: Neural_Net_residual() evaluates the neural network residual. Discussion: The interface is designed to correspond to the function minimizer ppc_nelder_mead(). Licensing: This code is distributed under the MIT license. Modified: 16 June 2024 Author: Original C version by Rouben Rostamian. This version by John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: double *weights: the weights. int nweights: the number of weights. void *params: associated parameters. Output: Neural_Net_residual: residual associated with this set of weights. */ { struct Neural_Net_PDE *nn = params; /* nu is the number of training points. */ int nu = nn->nu; double sum = 0.0; for ( int i = 0; i < nn->nweights; i++ ) { nn->weights[i] = weights[i]; } for ( int i = 0; i < nu; i++ ) { double x = nn->training_points[i][0]; double y = nn->training_points[i][1]; double r = residual_at_x_y ( nn, x, y ); sum = sum + r * r; } return sum; } /******************************************************************************/ double residual_at_x_y ( struct Neural_Net_PDE *nn, double x, double y ) /******************************************************************************/ /* Purpose: residual_at_x_y() evaluates the PDE residual at the point (x,y). Discussion: R(x,y) = F ( x, y, u, ux, uy, uxx, uxy, uyy ) Licensing: This code is distributed under the MIT license. Modified: 16 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 Input: struct Neural_Net_PDE *nn: describes the neural network. double x, y: the evaluation point. Output: double value: the PDE residual at x, y. */ { double u; double ux; double uy; double uxx; double uxy; double uyy; double value; nn->phi_func ( nn, x, y ); Neural_Net_eval ( nn, x, y ); u = nn->phi[0][0] * nn->N[0][0]; ux = nn->phi[1][0] * nn->N[0][0] + nn->phi[0][0] * nn->N[1][0]; uy = nn->phi[0][1] * nn->N[0][0] + nn->phi[0][0] * nn->N[0][1]; uxx = nn->phi[2][0] * nn->N[0][0] + 2.0 * nn->phi[1][0] * nn->N[1][0] + nn->phi[0][0] * nn->N[2][0]; uxy = nn->phi[1][1] * nn->N[0][0] + nn->phi[1][0] * nn->N[0][1] + nn->phi[0][1] * nn->N[1][0] + nn->phi[0][0] * nn->N[1][1]; uyy = nn->phi[0][2] * nn->N[0][0] + 2.0 * nn->phi[0][1] * nn->N[0][1] + nn->phi[0][0] * nn->N[0][2]; value = nn->PDE ( x, y, u, ux, uy, uxx, uxy, uyy ); return value; } /******************************************************************************/ void show_usage ( char *progname ) /******************************************************************************/ /* Purpose: show_usage() reminds the user how to call the program. Licensing: This code is distributed under the MIT license. Modified: 16 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 */ { printf ( "\n" ); printf ( "Usage:\n" ); printf ( " ./%s q nu\n", progname ); printf ( " q: number of units in the hidden layer (1 <= q)\n" ); printf ( " nu : number of training points (1 <= nu)\n" ); return; } /******************************************************************************/ void sigmoid ( double x, double *sigma ) /******************************************************************************/ /* Purpose: sigmoid() evaluates the sigmoidal function and its first two derivatives. Discussion: sigmoid(x) = 1 / ( 1 + exp(-x) ) Licensing: This code is distributed under the MIT license. Modified: 15 June 2024 Author: John Burkardt Reference: Rouben Rostamian, Programming Projects in C for Students of Engineering, Science, and Mathematics, SIAM, 2014, ISBN: 978-1-611973-49-5 */ { double value; value = 1.0 / ( 1.0 + exp ( - x ) ); sigma[0] = value; sigma[1] = value * ( 1.0 - value ); sigma[2] = value * ( 1.0 - value ) * ( 1.0 - 2.0 * value ); return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 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 }