# include # include # include # include "ppc_array.h" # include "ppc_nelder_mead.h" # include "ppc_neural_pde.h" # include "ppc_xmalloc.h" int main ( int argc, char **argv ); double square_hole_exact_sol ( double x, double y ); void square_hole_phi ( struct Neural_Net_PDE *nn, double x, double y ); double square_hole_pde ( double x, double y, double u, double u_x, double u_y, double u_xx, double u_xy, double u_yy ); /******************************************************************************/ int main ( int argc, char **argv ) /******************************************************************************/ /* Purpose: ppc_neural_pde_square_hole() tests ppc_neural_pde() on a square region with hole. Licensing: This code is distributed under the MIT license. Modified: 17 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 */ { char *endptr; double **training_points; timestamp ( ); printf ( "\n" ); printf ( "ppc_neural_pde_square_hole():\n" ); printf ( " ppc_neural_pde() models a partial differential equation\n" ); printf ( " in the form of a boundary value problem\n" ); printf ( " over a square region with a hole,\n" ); printf ( " solved using a neural network.\n" ); if ( argc != 3 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } /* Number of units in the hidden layer. */ int q = strtol ( argv[1], &endptr, 10 ); if ( *endptr != '\0' || q < 1 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } /* Number of training points. */ int nu = strtol ( argv[2], &endptr, 10 ); if ( *endptr != '\0' || nu < 1 ) { show_usage ( argv[0] ); return EXIT_FAILURE; } /* Create space for training points. */ make_matrix ( training_points, nu, 2 ); /* Define the neural network. */ struct Neural_Net_PDE nn = { .PDE = square_hole_pde, .phi_func = square_hole_phi, .bb_xrange = { -1.0, 1.0 }, .bb_yrange = { -1.0, 1.0 }, .q = q, .nu = nu, .training_points = training_points, .exact_sol = square_hole_exact_sol, .geomview_out = "square_hole.gv", .maple_out = "square_hole.m", .matlab_out = "square_hole.mpl", }; /* Generate nu training points. Do this by randomly picking points in the bounding box, and then checking that 0 < phi(x,y). */ int i = 0; while ( i < nu ) { double r = ( double ) rand() / RAND_MAX; double s = ( double ) rand() / RAND_MAX; double x = ( 1.0 - r ) * nn.bb_xrange[0] + r * nn.bb_xrange[1]; double y = ( 1.0 - s ) * nn.bb_yrange[0] + s * nn.bb_yrange[1]; nn.phi_func ( &nn, x, y ); if ( 0 < nn.phi[0][0] ) { nn.training_points[i][0] = x; nn.training_points[i][1] = y; i++; } } /* Initialize the neural net. */ Neural_Net_init ( &nn ); /* Create a structure for Nelder-Mead */ struct nelder_mead NM = { .f = Neural_Net_residual, .n = nn.nweights, .s = NULL, .x = nn.weights, .h = 0.1, .tol = 1e-5, .maxevals = 100000, .params = &nn, }; printf ( " weights before training:\n" ); print_vector ( "%7.3f ", nn.weights, nn.nweights ); /* Call nelder-mead to minimize the error function. */ int evalcount = nelder_mead ( &NM ); /* Report the results. */ if ( NM.maxevals < evalcount ) { printf ( " Nelder-Mead: No convergence after %d " "function evaluation\n", evalcount ); return EXIT_FAILURE; } printf ( " Nelder-Mead: Converged after %d function evaluations\n", evalcount ); printf ( " Nelder-Mead: Neural network’s residual error = %g\n", NM.minval ); printf ( " weights after training:\n"); print_vector ( "%7.3f ", nn.weights, nn.nweights ); if ( nn.exact_sol != NULL ) { printf ( " Error versus the ODE’s exact solution = %g\n", Neural_Net_error_vs_exact ( &nn, 50, 50 ) ); } Neural_Net_plot_with_maple ( &nn, 50, 50, "ppc_neural_pde_square_hole.mpl" ); Neural_Net_plot_with_matlab ( &nn, 50, 50, "ppc_neural_pde_square_hole.m" ); Neural_Net_end ( &nn ); free_vector ( training_points ); /* Terminate. */ printf ( "\n" ); printf ( "ppc_neural_ode_square_hole():\n" ); printf ( " Normal end of execution.\n" ); timestamp ( ); return EXIT_SUCCESS; } /******************************************************************************/ double square_hole_exact_sol ( double x, double y ) /******************************************************************************/ /* Purpose: square_hole_exact_sol() evaluates the exact solution. Licensing: This code is distributed under the MIT license. Discussion: For the exact solution we pick u(x,y) = (2+x) * (1-x^2) * (1-y^2) * (x^2+y^2-1/4) Modified: 17 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: double x, y: the evaluation point. Output: double value: the exact solution of the BVP, evaluated at (x,y). */ { double value; value = 2.0 * pow ( x, 5 ) + 4.0 * pow ( x, 4 ) + 0.5 * ( 64.0 * y * y - 49.0 ) * x * x * x + 0.5 * ( 96.0 * y * y - 66.0 ) * x * x + 0.5 * ( 12.0 * y * y * y - 51.0 * y * y + 20.0 ) * x + 4.0 * pow ( y, 4 ) - 33.0 * y * y + 10.0; return value; } /******************************************************************************/ void square_hole_phi ( struct Neural_Net_PDE *nn, double x, double y ) /******************************************************************************/ /* Purpose: square_hole_phi() indicates whether a point is inside the region. Discussion: We want to choose a function phi(x,y) which distinguishes between points (x,y) in the bounding box. phi(x,y) < 0 means these points are not part of the region. The cut-off function is phi(x,y) = (1-x^2) * (1-y^2) * ( x^2 + y^2 - 1/4 ). The derivatives of phi where calculated symbolically in Maple. The result was translated to C code also in Maple. Licensing: This code is distributed under the MIT license. Modified: 17 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: a structure defining a neural network. double x, y: the evaluation point. Output: struct Neural_Net_PDE *nn: the structure has been updated with phi[0][0], phi[1][0], phi[0][1], phi[2][0], phi[1][1], phi[0][2]. */ { double t1 = x * x; double t2 = -t1 + 0.1e1; double t4 = y * y; double t5 = -t4 + 0.1e1; double t6 = t1 + t4 - 0.1e1 / 0.4e1; double t7 = t5 * t6; double t9 = x * t5; double t11 = t2 * t5; double t15 = t2 * y; double t23 = 0.2e1 * t11; nn->phi[0][0] = t2 * t7; nn->phi[1][0] = 0.2e1 * (t11 * x - t9 * t6); nn->phi[0][1] = 0.2e1 * (t11 * y - t15 * t6); nn->phi[2][0] = -0.8e1 * t1 * t5 + t23 - 0.2e1 * t7; nn->phi[1][1] = 0.4e1 * (x * y * t6 - t15 * x - t9 * y); nn->phi[0][2] = -0.8e1 * t2 * t4 - 0.2e1 * t2 * t6 + t23; return; } /******************************************************************************/ double square_hole_pde ( double x, double y, double u, double u_x, double u_y, double u_xx, double u_xy, double u_yy ) /******************************************************************************/ /* Purpose: square_hole_pde() evaluates the PDE residual for the square problem. Discussion: For the exact solution we pick u(x,y) = (2+x) * (1-x^2) * (1-y^2) * (x^2+y^2-1/4) and then we "reverse engineer" the PDE by picking f = u_xx + u_yy. The calculation of f was done symbolically in Maple. Licensing: This code is distributed under the MIT license. Modified: 17 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: double x, y: the evaluation point. Output: double value: the value of the PDE residual, uxx + uyy - f. */ { double f; double t1; double t17; double t2; double t22; double t6; double value; t1 = x * x; t2 = t1 * t1; t6 = y * y; t17 = sqrt ( t1 + t6 ); t22 = t6 * t6; f = 0.2e1 / t17 * (t17 * (t2 + 0.3e1 * y * t1 * x + t1 * (0.6e1 * t6 - 0.7e1) + x * (0.3e1 * t6 * y - 0.6e1 * y) - t6 + 0.1e1) + x * (t2 + t1 * (0.19e2 / 0.2e1 * t6 - 0.15e2 / 0.2e1) + 0.3e1 * t22 - 0.15e2 / 0.2e1 * t6 + 0.3e1 / 0.2e1)); value = u_xx + u_yy - f; return value; }