# include # include # include "minimal_surface_exact.h" /******************************************************************************/ void minimal_surface_catenoid_exact ( int n, double X[], double Y[], double a, double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], double Uyy[] ) /******************************************************************************/ /* Purpose: minimal_surface_catenoid_exact() evaluates catenoid minimal surface solution U(X,Y). Discussion: The minimal surface equation describes a function with zero mean curvature. The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The helicoid solution is: U(X,Y) = acosh ( a sqrt(X^2+Y^2) ) / a Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Reference: John D Cook, Closed-form minimal surface solutions, https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/ Posted 31 March 2025. Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a: a shape parameter. Output: double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): the solution, first and second partial derivatives evaluated at (X,Y). */ { int i; double n1; double n2; for ( i = 0; i < n; i++ ) { n1 = X[i] * X[i] + Y[i] * Y[i]; n2 = sqrt ( n1 ); U[i] = acosh ( a * n2 ) / a; Ux[i] = X[i] / ( n2 * sqrt ( a * a * n1 - 1.0 ) ); Uy[i] = Y[i] / ( n2 * sqrt ( a * a * n1 - 1.0 ) ); Uxx[i] = - ( a * a * X[i] * X[i] / ( a * a * n1 - 1.0 ) + X[i] * X[i] / n1 - 1.0 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) ); Uxy[i] = - X[i] * Y[i] * ( a * a / ( a * a * n1 - 1.0 ) + 1.0 / n1 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) ); Uyy[i] = - ( a * a * Y[i] * Y[i] / ( a * a * n1 - 1.0 ) + Y[i] * Y[i] / n1 - 1.0 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) ); } return; } /******************************************************************************/ double *minimal_surface_catenoid_residual ( int n, double X[], double Y[], double a ) /******************************************************************************/ /* Purpose: minimal_surface_catenoid_residual(): minimal surface residual for catenoid solution. Discussion: The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The catenoid solution is: U(X,Y) = acosh ( a * sqrt ( X^2 + Y^2 ) ) / a Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a: a shape parameter. Output: double R(n): the residual evaluated at the points (X,Y). */ { int i; double *R; double *U; double *Ux; double *Uxx; double *Uxy; double *Uy; double *Uyy; U = ( double * ) malloc ( n * sizeof ( double ) ); Ux = ( double * ) malloc ( n * sizeof ( double ) ); Uy = ( double * ) malloc ( n * sizeof ( double ) ); Uxx = ( double * ) malloc ( n * sizeof ( double ) ); Uxy = ( double * ) malloc ( n * sizeof ( double ) ); Uyy = ( double * ) malloc ( n * sizeof ( double ) ); minimal_surface_catenoid_exact ( n, X, Y, a, U, Ux, Uy, Uxx, Uxy, Uyy ); R = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] - 2.0 * Ux[i] * Uxy[i] * Uy[i] + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i]; } free ( U ); free ( Ux ); free ( Uy ); free ( Uxx ); free ( Uxy ); free ( Uyy ); return R; } /******************************************************************************/ void minimal_surface_helicoid_exact ( int n, double X[], double Y[], double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], double Uyy[] ) /******************************************************************************/ /* Purpose: minimal_surface_helicoid_exact() evaluates helicoid minimal surface solution U(X,Y). Discussion: The minimal surface equation describes a function with zero mean curvature. The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The helicoid solution is: U(X,Y) = arctan ( X / Y ) Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Reference: John D Cook, Closed-form minimal surface solutions, https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/ Posted 31 March 2025. Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. Output: double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): the solution, first and second partial derivatives evaluated at (X,Y). */ { int i; for ( i = 0; i < n; i++ ) { U[i] = atan ( X[i] / Y[i] ); Ux[i] = 1.0 / ( Y[i] * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ); Uy[i] = - X[i] / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ); Uxx[i] = - 2.0 * X[i] / ( pow ( Y[i], 3 ) * pow ( ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ), 2 ) ); Uxy[i] = ( 2.0 * pow ( X[i], 2 ) / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ) - 1.0 ) / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ); Uyy[i] = 2.0 * X[i] * ( - pow ( X[i], 2 ) / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ) + 1.0 ) / ( pow ( Y[i], 3 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ); } return; } /******************************************************************************/ double *minimal_surface_helicoid_residual ( int n, double X[], double Y[] ) /******************************************************************************/ /* Purpose: minimal_surface_helicoid_residual(): minimal surface residual for helicoid solution. Discussion: The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The helicoid solution is: U(X,Y) = arctan ( X / Y ) Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. Output: double R(n): the residual evaluated at the points (X,Y). */ { int i; double *R; double *U; double *Ux; double *Uxx; double *Uxy; double *Uy; double *Uyy; U = ( double * ) malloc ( n * sizeof ( double ) ); Ux = ( double * ) malloc ( n * sizeof ( double ) ); Uy = ( double * ) malloc ( n * sizeof ( double ) ); Uxx = ( double * ) malloc ( n * sizeof ( double ) ); Uxy = ( double * ) malloc ( n * sizeof ( double ) ); Uyy = ( double * ) malloc ( n * sizeof ( double ) ); minimal_surface_helicoid_exact ( n, X, Y, U, Ux, Uy, Uxx, Uxy, Uyy ); R = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] - 2.0 * Ux[i] * Uxy[i] * Uy[i] + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i]; } free ( U ); free ( Ux ); free ( Uy ); free ( Uxx ); free ( Uxy ); free ( Uyy ); return R; } /******************************************************************************/ void minimal_surface_linear_exact ( int n, double X[], double Y[], double a, double b, double c, double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], double Uyy[] ) /******************************************************************************/ /* Purpose: minimal_surface_linear_exact() evaluates linear minimal surface solution U(X,Y). Discussion: The minimal surface equation describes a function with zero mean curvature. The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The linear solution is: U(X,Y) = a * X + b * Y + C Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Reference: John D Cook, Closed-form minimal surface solutions, https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/ Posted 31 March 2025. Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a, b, c: parameters. Output: double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): the solution, first and second partial derivatives evaluated at (X,Y). */ { int i; for ( i = 0; i < n; i++ ) { U[i] = a * X[i] + b * Y[i] + c; Ux[i] = a; Uy[i] = b; Uxx[i] = 0.0; Uxy[i] = 0.0; Uyy[i] = 0.0; } return; } /******************************************************************************/ double *minimal_surface_linear_residual ( int n, double X[], double Y[], double a, double b, double c ) /******************************************************************************/ /* Purpose: minimal_surface_linear_residual(): minimal surface residual for linear solution. Discussion: The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The linear solution is: U(X,Y) = a * X + b * Y + C Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a, b, c: parameters. Output: double R(n): the residual evaluated at the points (X,Y). */ { int i; double *R; double *U; double *Ux; double *Uxx; double *Uxy; double *Uy; double *Uyy; U = ( double * ) malloc ( n * sizeof ( double ) ); Ux = ( double * ) malloc ( n * sizeof ( double ) ); Uy = ( double * ) malloc ( n * sizeof ( double ) ); Uxx = ( double * ) malloc ( n * sizeof ( double ) ); Uxy = ( double * ) malloc ( n * sizeof ( double ) ); Uyy = ( double * ) malloc ( n * sizeof ( double ) ); minimal_surface_linear_exact ( n, X, Y, a, b, c, U, Ux, Uy, Uxx, Uxy, Uyy ); R = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] - 2.0 * Ux[i] * Uxy[i] * Uy[i] + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i]; } free ( U ); free ( Ux ); free ( Uy ); free ( Uxx ); free ( Uxy ); free ( Uyy ); return R; } /******************************************************************************/ void minimal_surface_scherk_exact ( int n, double X[], double Y[], double a, double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], double Uyy[] ) /******************************************************************************/ /* Purpose: minimal_surface_scherk_exact() evaluates Scherk minimal surface solution U(X,Y). Discussion: The minimal surface equation describes a function with zero mean curvature. The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The Scherk solution is: U(X,Y) = log ( cos ( a Y ) / cos ( a X ) ) / a Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Reference: John D Cook, Closed-form minimal surface solutions, https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/ Posted 31 March 2025. Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a: a shape parameter. Output: double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): the solution, first and second partial derivatives evaluated at (X,Y). */ { int i; for ( i = 0; i < n; i++ ) { U[i] = log ( cos ( a * Y[i] ) / cos ( a * X[i] ) ) / a; Ux[i] = sin ( a * X[i] ) / cos ( a * X[i] ); Uy[i] = - sin ( a * Y[i] ) / cos ( a * Y[i] ); Uxx[i] = a * ( pow ( sin ( a * X[i] ), 2 ) / pow ( cos ( a * X[i] ), 2 ) + 1.0 ); Uxy[i] = 0.0; Uyy[i] = - a * ( pow ( sin ( a * Y[i] ), 2 ) / pow ( cos ( a * Y[i] ), 2 ) + 1.0 ); } return; } /******************************************************************************/ double *minimal_surface_scherk_residual ( int n, double X[], double Y[], double a ) /******************************************************************************/ /* Purpose: minimal_surface_scherk_residual(): minimal surface residual for Scherk solution. Discussion: The equation is: (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0 The Scherk solution is: U(X,Y) = log ( cos ( a * Y ) / cos ( a * X ) ) / a Licensing: This code is distributed under the MIT license. Modified: 02 June 2025 Author: John Burkardt Input: int n: the number of points. double X(n), Y(n): the coordinates of the points. double a: a shape parameter. Output: double R(n): the residual evaluated at the points (X,Y). */ { int i; double *R; double *U; double *Ux; double *Uxx; double *Uxy; double *Uy; double *Uyy; U = ( double * ) malloc ( n * sizeof ( double ) ); Ux = ( double * ) malloc ( n * sizeof ( double ) ); Uy = ( double * ) malloc ( n * sizeof ( double ) ); Uxx = ( double * ) malloc ( n * sizeof ( double ) ); Uxy = ( double * ) malloc ( n * sizeof ( double ) ); Uyy = ( double * ) malloc ( n * sizeof ( double ) ); minimal_surface_scherk_exact ( n, X, Y, a, U, Ux, Uy, Uxx, Uxy, Uyy ); R = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] - 2.0 * Ux[i] * Uxy[i] * Uy[i] + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i]; } free ( U ); free ( Ux ); free ( Uy ); free ( Uxx ); free ( Uxy ); free ( Uyy ); return R; }