# 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 ell_exact_sol ( double x, double y ); void ell_phi ( struct Neural_Net_PDE *nn, double x, double y ); double ell_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_ell() tests ppc_neural_pde() on the L-shaped region. Discussion: Seek a solution in an L-shaped region of the boundary value problem: d2 u/dx2 + d2 u/dy2 = f(x,y) where u(x,y) = 0 on the boundary and f(x,y) = 32 ( x^2 + y^2 - x - y ) for which the exact solution is u(x,y) = 16 * x * ( x - 1 ) * y * ( y - 1 ) 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_ell():\n" ); printf ( " ppc_neural_pde() models a partial differential equation\n" ); printf ( " in the form of a boundary value problem over an L-shaped region\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 = ell_pde, .phi_func = ell_phi, .bb_xrange = { -1.0, 1.0 }, .bb_yrange = { -1.0, 1.0 }, .q = q, .nu = nu, .training_points = training_points, .exact_sol = ell_exact_sol, .geomview_out = "ell.gv", .maple_out = "ell.m", .matlab_out = "ell.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_ell.mpl" ); Neural_Net_plot_with_matlab ( &nn, 50, 50, "ppc_neural_pde_ell.m" ); Neural_Net_end ( &nn ); free_vector ( training_points ); /* Terminate. */ printf ( "\n" ); printf ( "ppc_neural_ode_ell():\n" ); printf ( " Normal end of execution.\n" ); timestamp ( ); return EXIT_SUCCESS; } /******************************************************************************/ double ell_exact_sol ( double x, double y ) /******************************************************************************/ /* Purpose: ell_exact_sol() evaluates the exact solution. 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 exact solution of the BVP, evaluated at (x,y). */ { double value; value = x * pow ( 1.0 - x, 2 ) * pow ( 1.0 - y, 2 ) * ( x + y + sqrt ( x * x + y * y ) ); return value; } /******************************************************************************/ void ell_phi ( struct Neural_Net_PDE *nn, double x, double y ) /******************************************************************************/ /* Purpose: ell_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 for the L-shaped region is phi(x,y) = (1-x^2)*(1-y^2)*(x + y + sqrt(x^2+y^2)). 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; double t7 = sqrt(t6); double t8 = t7 + y + x; double t9 = t8 * t5; double t11 = t5 * x; double t14 = t5 * t2; double t15 = 0.1e1 / t7; double t17 = t15 * x + 0.1e1; double t21 = y * t2; double t25 = t15 * y + 0.1e1; double t33 = 0.1e1 / t7 / t6; nn->phi[0][0] = t9 * t2; nn->phi[1][0] = 0.2e1 * t8 * t11 + t17 * t14; nn->phi[0][1] = t25 * t14 + 0.2e1 * t8 * t21; nn->phi[2][0] = 0.2e1 * t9 + 0.4e1 * t17 * t11 + (-t33 * t1 + t15) * t14; nn->phi[1][1] = -y * t33 * x * t14 + 0.4e1 * t8 * y * x + 0.2e1 * t25 * t11 + 0.2e1 * t17 * t21; nn->phi[0][2] = 0.2e1 * t8 * t2 + 0.4e1 * t25 * t21 + (-t33 * t4 + t15) * t14; return; } /******************************************************************************/ double ell_pde ( double x, double y, double u, double u_x, double u_y, double u_xx, double u_xy, double u_yy ) /******************************************************************************/ /* Purpose: ell_pde() evaluates the PDE residual for the L-shaped region problem. Discussion: For the exact solution we pick u(x,t) = x*phi(x,y), 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 t1 = x * x; double t2 = t1 * t1; double t6 = y * y; double t17 = sqrt ( t1 + t6 ); double t22 = t6 * t6; double 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)); double value; value = u_xx + u_yy - f; return value; }