# include # include # include # include # include "fem1d_bvp_linear.h" int main ( ); void test00 ( ); double a00 ( double x ); double c00 ( double x ); double exact00 ( double x ); double exact_ux00 ( double x ); double f00 ( double x ); void test01 ( ); double a1 ( double x ); double c1 ( double x ); double exact1 ( double x ); double exact_ux1 ( double x ); double f1 ( double x ); void test02 ( ); double a2 ( double x ); double c2 ( double x ); double exact2 ( double x ); double exact_ux2 ( double x ); double f2 ( double x ); void test03 ( ); double a3 ( double x ); double c3 ( double x ); double exact3 ( double x ); double exact_ux3 ( double x ); double f3 ( double x ); void test04 ( ); double a4 ( double x ); double c4 ( double x ); double exact4 ( double x ); double exact_ux4 ( double x ); double f4 ( double x ); void test05 ( ); double a5 ( double x ); double c5 ( double x ); double exact5 ( double x ); double exact_ux5 ( double x ); double f5 ( double x ); void test06 ( ); double a6 ( double x ); double c6 ( double x ); double exact6 ( double x ); double exact_ux6 ( double x ); double f6 ( double x ); void test07 ( ); double a7 ( double x ); double c7 ( double x ); double exact7 ( double x ); double exact_ux7 ( double x ); double f7 ( double x ); void test08 ( ); double a8 ( double x ); double c8 ( double x ); double exact8 ( double x ); double exact_ux8 ( double x ); double f8 ( double x ); void test09 ( ); double a9 ( double x ); double c9 ( double x ); double exact9 ( double x ); double exact_ux9 ( double x ); double f9 ( double x ); void test10 ( ); double a10 ( double x ); double c10 ( double x ); double exact10 ( double x ); double exact_ux10 ( double x ); double f10 ( double x ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: fem1d_bvp_linear_test() tests fem1d_bvp_linear(). Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "FEM1D_BVP_LINEAR_TEST\n" ); printf ( " C version\n" ); printf ( " Test the FEM1D_BVP_LINEAR library.\n" ); test00 ( ); test01 ( ); test02 ( ); test03 ( ); test04 ( ); test05 ( ); test06 ( ); test07 ( ); test08 ( ); test09 ( ); test10 ( ); /* Terminate. */ printf ( "\n" ); printf ( "FEM1D_BVP_LINEAR_TEST\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void test00 ( ) /******************************************************************************/ /* Purpose: TEST00 carries out test case #0. Discussion: - uxx + u = x for 0 < x < 1 u(0) = u(1) = 0 exact = x - sinh(x) / sinh(1) exact' = 1 - cosh(x) / sinh(1) Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST00\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A(X) = 1.0\n" ); printf ( " C(X) = 1.0\n" ); printf ( " F(X) = X\n" ); printf ( " U(X) = X - SINH(X) / SINH(1)\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a00, c00, f00, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact00 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact00 ); e2 = l2_error_linear ( n, x, u, exact00 ); h1s = h1s_error_linear ( n, x, u, exact_ux00 ); mx = max_error_linear ( n, x, u, exact00 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a00 ( double x ) /******************************************************************************/ /* Purpose: A00 evaluates A function #0. Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A00, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c00 ( double x ) /******************************************************************************/ /* Purpose: C00 evaluates C function #0. Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C00, the value of C(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double exact00 ( double x ) /******************************************************************************/ /* Purpose: EXACT00 evaluates exact solution #0. Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT00, the value of U(X). */ { double value; value = x - sinh ( x ) / sinh ( 1.0 ); return value; } /******************************************************************************/ double exact_ux00 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX00 evaluates the derivative of exact solution #0. Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX00, the value of dUdX(X). */ { double value; value = 1.0 - cosh ( x ) / sinh ( 1.0 ); return value; } /******************************************************************************/ double f00 ( double x ) /******************************************************************************/ /* Purpose: F00 evaluates right hand side function #0. Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F00, the value of F(X). */ { double value; value = x; return value; } /******************************************************************************/ void test01 ( ) /******************************************************************************/ /* Purpose: TEST01 carries out test case #1. Discussion: Use A1, C1, F1, EXACT1, EXACT_UX1. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST01\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A1(X) = 1.0\n" ); printf ( " C1(X) = 0.0\n" ); printf ( " F1(X) = X * ( X + 3 ) * exp ( X )\n" ); printf ( " U1(X) = X * ( 1 - X ) * exp ( X )\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a1, c1, f1, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact1 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact1 ); e2 = l2_error_linear ( n, x, u, exact1 ); h1s = h1s_error_linear ( n, x, u, exact_ux1 ); mx = max_error_linear ( n, x, u, exact1 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a1 ( double x ) /******************************************************************************/ /* Purpose: A1 evaluates A function #1. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A1, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c1 ( double x ) /******************************************************************************/ /* Purpose: C1 evaluates C function #1. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C1, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact1 ( double x ) /******************************************************************************/ /* Purpose: EXACT1 evaluates exact solution #1. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT1, the value of U(X). */ { double value; value = x * ( 1.0 - x ) * exp ( x ); return value; } /******************************************************************************/ double exact_ux1 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX1 evaluates the derivative of exact solution #1. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX1, the value of dUdX(X). */ { double value; value = ( 1.0 - x - x * x ) * exp ( x ); return value; } /******************************************************************************/ double f1 ( double x ) /******************************************************************************/ /* Purpose: F1 evaluates right hand side function #1. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F1, the value of F(X). */ { double value; value = x * ( x + 3.0 ) * exp ( x ); return value; } /******************************************************************************/ void test02 ( ) /******************************************************************************/ /* Purpose: TEST02 carries out test case #2. Discussion: Use A2, C2, F2, EXACT2, EXACT_UX2. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST02\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A2(X) = 1.0\n" ); printf ( " C2(X) = 2.0\n" ); printf ( " F2(X) = X * ( 5 - X ) * exp ( X )\n" ); printf ( " U2(X) = X * ( 1 - X ) * exp ( X )\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a2, c2, f2, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact2 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact2 ); e2 = l2_error_linear ( n, x, u, exact2 ); h1s = h1s_error_linear ( n, x, u, exact_ux2 ); mx = max_error_linear ( n, x, u, exact2 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a2 ( double x ) /******************************************************************************/ /* Purpose: A2 evaluates A function #2. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A2, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c2 ( double x ) /******************************************************************************/ /* Purpose: C2 evaluates C function #2. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C2, the value of C(X). */ { double value; value = 2.0; return value; } /******************************************************************************/ double exact2 ( double x ) /******************************************************************************/ /* Purpose: EXACT2 evaluates exact solution #2. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT2, the value of U(X). */ { double value; value = x * ( 1.0 - x ) * exp ( x ); return value; } /******************************************************************************/ double exact_ux2 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX2 evaluates the derivative of exact solution #2. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX2, the value of dUdX(X). */ { double value; value = ( 1.0 - x - x * x ) * exp ( x ); return value; } /******************************************************************************/ double f2 ( double x ) /******************************************************************************/ /* Purpose: F2 evaluates right hand side function #2. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F2, the value of F(X). */ { double value; value = x * ( 5.0 - x ) * exp ( x ); return value; } /******************************************************************************/ void test03 ( ) /******************************************************************************/ /* Purpose: TEST03 carries out test case #3. Discussion: Use A3, C3, F3, EXACT3, EXACT_UX3. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST03\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A3(X) = 1.0\n" ); printf ( " C3(X) = 2.0 * X\n" ); printf ( " F3(X) = - X * ( 2 * X * X - 3 * X - 3 ) * exp ( X )\n" ); printf ( " U3(X) = X * ( 1 - X ) * exp ( X )\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a3, c3, f3, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact3 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact3 ); e2 = l2_error_linear ( n, x, u, exact3 ); h1s = h1s_error_linear ( n, x, u, exact_ux3 ); mx = max_error_linear ( n, x, u, exact3 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a3 ( double x ) /******************************************************************************/ /* Purpose: A3 evaluates A function #3. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A3, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c3 ( double x ) /******************************************************************************/ /* Purpose: C3 evaluates C function #3. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C3, the value of C(X). */ { double value; value = 2.0 * x; return value; } /******************************************************************************/ double exact3 ( double x ) /******************************************************************************/ /* Purpose: EXACT3 evaluates exact solution #3. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT3, the value of U(X). */ { double value; value = x * ( 1.0 - x ) * exp ( x ); return value; } /******************************************************************************/ double exact_ux3 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX3 evaluates the derivative of exact solution #3. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX3, the value of dUdX(X). */ { double value; value = ( 1.0 - x - x * x ) * exp ( x ); return value; } /******************************************************************************/ double f3 ( double x ) /******************************************************************************/ /* Purpose: F3 evaluates right hand side function #3. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F3, the value of F(X). */ { double value; value = - x * ( 2.0 * x * x - 3.0 * x - 3.0 ) * exp ( x ); return value; } /******************************************************************************/ void test04 ( ) /******************************************************************************/ /* Purpose: TEST04 carries out test case #4. Discussion: Use A4, C4, F4, EXACT4, EXACT_UX4. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST04\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A4(X) = 1.0 + X * X\n" ); printf ( " C4(X) = 0.0\n" ); printf ( " F4(X) = ( X + 3 X^2 + 5 X^3 + X^4 ) * exp ( X )\n" ); printf ( " U4(X) = X * ( 1 - X ) * exp ( X )\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a4, c4, f4, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact4 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact4 ); e2 = l2_error_linear ( n, x, u, exact4 ); h1s = h1s_error_linear ( n, x, u, exact_ux4 ); mx = max_error_linear ( n, x, u, exact4 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a4 ( double x ) /******************************************************************************/ /* Purpose: A4 evaluates A function #4. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A4, the value of A(X). */ { double value; value = 1.0 + x * x; return value; } /******************************************************************************/ double c4 ( double x ) /******************************************************************************/ /* Purpose: C4 evaluates C function #4. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C4, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact4 ( double x ) /******************************************************************************/ /* Purpose: EXACT4 evaluates exact solution #4. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT4, the value of U(X). */ { double value; value = x * ( 1.0 - x ) * exp ( x ); return value; } /******************************************************************************/ double exact_ux4 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX4 evaluates the derivative of exact solution #4. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX4, the value of dUdX(X). */ { double value; value = ( 1.0 - x - x * x ) * exp ( x ); return value; } /******************************************************************************/ double f4 ( double x ) /******************************************************************************/ /* Purpose: F4 evaluates right hand side function #4. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F4, the value of F(X). */ { double value; value = ( x + 3.0 * x * x + 5.0 * x * x * x + x * x * x * x ) * exp ( x ); return value; } /******************************************************************************/ void test05 ( ) /******************************************************************************/ /* Purpose: TEST05 carries out test case #5. Discussion: Use A5, C5, F5, EXACT5, EXACT_UX5. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST05\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A5(X) = 1.0 + X * X for X <= 1/3\n" ); printf ( " = 7/9 + X for 1/3 < X\n" ); printf ( " C5(X) = 0.0\n" ); printf ( " F5(X) = ( X + 3 X^2 + 5 X^3 + X^4 ) * exp ( X )\n" ); printf ( " for X <= 1/3\n" ); printf ( " = ( - 1 + 10/3 X + 43/9 X^2 + X^3 ) .* exp ( X )\n" ); printf ( " for 1/3 <= X\n" ); printf ( " U5(X) = X * ( 1 - X ) * exp ( X )\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a5, c5, f5, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact5 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact5 ); e2 = l2_error_linear ( n, x, u, exact5 ); h1s = h1s_error_linear ( n, x, u, exact_ux5 ); mx = max_error_linear ( n, x, u, exact5 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a5 ( double x ) /******************************************************************************/ /* Purpose: A5 evaluates A function #5. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A5, the value of A(X). */ { double value; if ( x <= 1.0 / 3.0 ) { value = 1.0 + x * x; } else { value = x + 7.0 / 9.0; } return value; } /******************************************************************************/ double c5 ( double x ) /******************************************************************************/ /* Purpose: C5 evaluates C function #5. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C5, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact5 ( double x ) /******************************************************************************/ /* Purpose: EXACT5 evaluates exact solution #5. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT5, the value of U(X). */ { double value; value = x * ( 1.0 - x ) * exp ( x ); return value; } /******************************************************************************/ double exact_ux5 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX5 evaluates the derivative of exact solution #5. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX5, the value of dUdX(X). */ { double value; value = ( 1.0 - x - x * x ) * exp ( x ); return value; } /******************************************************************************/ double f5 ( double x ) /******************************************************************************/ /* Purpose: F5 evaluates right hand side function #5. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F5, the value of F(X). */ { double value; if ( x <= 1.0 / 3.0 ) { value = ( x + 3.0 * x * x + 5.0 * x * x * x + x * x * x * x ) * exp ( x ); } else { value = ( - 1.0 + ( 10.0 / 3.0 ) * x + ( 43.0 / 9.0 ) * x * x + x * x * x ) * exp ( x ); } return value; } /******************************************************************************/ void test06 ( ) /******************************************************************************/ /* Purpose: TEST06 does an error analysis. Discussion: Use A6, C6, F6, EXACT6, EXACT_UX6. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt */ { int i; int n; double e1; double e2; double h1s; double mx; double *u; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST06\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A6(X) = 1.0\n" ); printf ( " C6(X) = 0.0\n" ); printf ( " F6(X) = pi*pi*sin(pi*X)\n" ); printf ( " U6(X) = sin(pi*x)\n" ); printf ( "\n" ); printf ( " Compute l1norm, L2norm and seminorm of error for various N.\n" ); printf ( "\n" ); printf ( " N l1 error L2 error Seminorm error Maxnorm error\n" ); printf ( "\n" ); n = 11; for ( i = 0; i <= 4; i++ ) { /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a6, c6, f6, x ); e1 = l1_error ( n, x, u, exact6 ); e2 = l2_error_linear ( n, x, u, exact6 ); h1s = h1s_error_linear ( n, x, u, exact_ux6 ); mx = max_error_linear ( n, x, u, exact6 ); printf ( " %4d %14g %14g %14g %14g\n", n, e1, e2, h1s, mx ); n = 2 * ( n - 1 ) + 1; free ( u ); free ( x ); } return; } /******************************************************************************/ double a6 ( double x ) /******************************************************************************/ /* Purpose: A6 evaluates A function #6. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A6, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c6 ( double x ) /******************************************************************************/ /* Purpose: C6 evaluates C function #6. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C6, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact6 ( double x ) /******************************************************************************/ /* Purpose: EXACT6 returns exact solution #6. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT6, the value of U(X). */ { const double pi = 3.141592653589793; double value; value = sin ( pi * x ); return value; } /******************************************************************************/ double exact_ux6 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX6 returns the derivative of exact solution #6. Licensing: This code is distributed under the MIT license. Modified: 14 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX6, the value of U(X). */ { const double pi = 3.141592653589793; double value; value = pi * cos ( pi * x ); return value; } /******************************************************************************/ double f6 ( double x ) /******************************************************************************/ /* Purpose: F6 evaluates right hand side function #6. Licensing: This code is distributed under the MIT license. Modified: 19 February 2012 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F6, the value of F(X). */ { static double pi = 3.141592653589793; double value; value = pi * pi * sin ( pi * x ); return value; } /******************************************************************************/ void test07 ( ) /******************************************************************************/ /* Purpose: TEST07 does an error analysis. Discussion: Use A7, C7, F7, EXACT7, EXACT_UX7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Reference: Eric Becker, Graham Carey, John Oden, Finite Elements, An Introduction, Volume I, Prentice-Hall, 1981, page 123-124, ISBN: 0133170578, LC: TA347.F5.B4. */ { int i; int n; double e1; double e2; double h1s; double mx; double *u; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST07\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " Becker/Carey/Oden Example\n" ); printf ( "\n" ); printf ( " Compute l1 norm, L2 norm and seminorm of error for various N.\n" ); printf ( "\n" ); printf ( " N l1 error L2 error Seminorm error Maxnorm error\n" ); printf ( "\n" ); n = 11; for ( i = 0; i <= 4; i++ ) { /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a7, c7, f7, x ); e1 = l1_error ( n, x, u, exact7 ); e2 = l2_error_linear ( n, x, u, exact7 ); h1s = h1s_error_linear ( n, x, u, exact_ux7 ); mx = max_error_linear ( n, x, u, exact7 ); printf ( " %4d %14g %14g %14g %14g\n", n, e1, e2, h1s, mx ); n = 2 * ( n - 1 ) + 1; free ( u ); free ( x ); } return; } /******************************************************************************/ double a7 ( double x ) /******************************************************************************/ /* Purpose: A7 evaluates A function #7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A7, the value of A(X). */ { double alpha; double value; double x0; alpha = 30.0; x0 = 1.0 / 3.0; value = 1.0 / alpha + alpha * pow ( x - x0, 2 ); return value; } /******************************************************************************/ double c7 ( double x ) /******************************************************************************/ /* Purpose: C7 evaluates C function #7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C7, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact7 ( double x ) /******************************************************************************/ /* Purpose: EXACT7 returns exact solution #7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT7, the value of U(X). */ { double alpha; double value; double x0; alpha = 30.0; x0 = 1.0 / 3.0; value = ( 1.0 - x ) * ( atan ( alpha * ( x - x0 ) ) + atan ( alpha * x0 ) ); return value; } /******************************************************************************/ double exact_ux7 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX7 returns the derivative of exact solution #7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX7, the value of U(X). */ { double alpha; double value; double x0; alpha = 30.0; x0 = 1.0 / 3.0; value = - atan ( alpha * ( x - x0 ) ) - atan ( alpha * x0 ) + ( 1.0 - x ) * alpha / ( 1.0 + alpha * alpha * pow ( x - x0, 2 ) ); return value; } /******************************************************************************/ double f7 ( double x ) /******************************************************************************/ /* Purpose: F7 evaluates right hand side function #7. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F7, the value of F(X). */ { double alpha; double value; double x0; alpha = 30.0; x0 = 1.0 / 3.0; value = 2.0 * ( 1.0 + alpha * ( x - x0 ) * ( atan ( alpha * ( x - x0 ) ) + atan ( alpha * x0 ) ) ); return value; } /******************************************************************************/ void test08 ( ) /******************************************************************************/ /* Purpose: TEST08 carries out test case #8. Discussion: Use A8, C8, F8, EXACT8, EXACT_UX8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST08\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A8(X) = 1.0\n" ); printf ( " C8(X) = 0.0\n" ); printf ( " F8(X) = X * ( X + 3 ) * exp ( X ), X <= 2/3\n" ); printf ( " = 2 * exp ( 2/3), 2/3 < X\n" ); printf ( " U8(X) = X * ( 1 - X ) * exp ( X ), X <= 2/3\n" ); printf ( " = X * ( 1 - X ) * exp ( 2/3 ), 2/3 < X\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a8, c8, f8, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact8 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact8 ); e2 = l2_error_linear ( n, x, u, exact8 ); h1s = h1s_error_linear ( n, x, u, exact_ux8 ); mx = max_error_linear ( n, x, u, exact8 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a8 ( double x ) /******************************************************************************/ /* Purpose: A8 evaluates A function #8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A8, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c8 ( double x ) /******************************************************************************/ /* Purpose: C8 evaluates C function #8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C8, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact8 ( double x ) /******************************************************************************/ /* Purpose: EXACT8 evaluates exact solution #8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT8, the value of U(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = x * ( 1.0 - x ) * exp ( x ); } else { value = x * ( 1.0 - x ) * exp ( 2.0 / 3.0 ); } return value; } /******************************************************************************/ double exact_ux8 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX8 evaluates the derivative of exact solution #8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX8, the value of dUdX(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = ( 1.0 - x - x * x ) * exp ( x ); } else { value = ( 1.0 - 2.0 * x ) * exp ( 2.0 / 3.0 ); } return value; } /******************************************************************************/ double f8 ( double x ) /******************************************************************************/ /* Purpose: F8 evaluates right hand side function #8. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F8, the value of F(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = x * ( x + 3.0 ) * exp ( x ); } else { value = 2.0 * exp ( 2.0 / 3.0 ); } return value; } /******************************************************************************/ void test09 ( ) /******************************************************************************/ /* Purpose: TEST09 carries out test case #9. Discussion: Use A9, C9, F9, EXACT9, EXACT_UX9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { int i; int n = 11; double e1; double e2; double h1s; double mx; double *u; double uexact; double *x; double x_first; double x_last; printf ( "\n" ); printf ( "TEST09\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A9(X) = 1.0\n" ); printf ( " C9(X) = 0.0\n" ); printf ( " F9(X) = X * ( X + 3 ) * exp ( X ), X <= 2/3\n" ); printf ( " = 2 * exp ( 2/3), 2/3 < X\n" ); printf ( " U9(X) = X * ( 1 - X ) * exp ( X ), X <= 2/3\n" ); printf ( " = X * ( 1 - X ), 2/3 < X\n" ); printf ( "\n" ); printf ( " Number of nodes = %d\n", n ); /* Geometry definitions. */ x_first = 0.0; x_last = 1.0; x = r8vec_linspace_new ( n, x_first, x_last ); u = fem1d_bvp_linear ( n, a9, c9, f9, x ); printf ( "\n" ); printf ( " I X U Uexact Error\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { uexact = exact9 ( x[i] ); printf ( " %4d %8f %14f %14f %14e\n", i, x[i], u[i], uexact, fabs ( u[i] - uexact ) ); } e1 = l1_error ( n, x, u, exact9 ); e2 = l2_error_linear ( n, x, u, exact9 ); h1s = h1s_error_linear ( n, x, u, exact_ux9 ); mx = max_error_linear ( n, x, u, exact9 ); printf ( "\n" ); printf ( " l1 norm of error = %g\n", e1 ); printf ( " L2 norm of error = %g\n", e2 ); printf ( " Seminorm of error = %g\n", h1s ); printf ( " Max norm of error = %g\n", mx ); free ( u ); free ( x ); return; } /******************************************************************************/ double a9 ( double x ) /******************************************************************************/ /* Purpose: A9 evaluates A function #9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A9, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c9 ( double x ) /******************************************************************************/ /* Purpose: C9 evaluates C function #9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C9, the value of C(X). */ { double value; value = 0.0; return value; } /******************************************************************************/ double exact9 ( double x ) /******************************************************************************/ /* Purpose: EXACT9 evaluates exact solution #9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT9, the value of U(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = x * ( 1.0 - x ) * exp ( x ); } else { value = x * ( 1.0 - x ); } return value; } /******************************************************************************/ double exact_ux9 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX9 evaluates the derivative of exact solution #9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX9, the value of dUdX(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = ( 1.0 - x - x * x ) * exp ( x ); } else { value = 1.0 - 2.0 * x; } return value; } /******************************************************************************/ double f9 ( double x ) /******************************************************************************/ /* Purpose: F9 evaluates right hand side function #9. Licensing: This code is distributed under the MIT license. Modified: 16 June 2014 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F9, the value of F(X). */ { double value; if ( x <= 2.0 / 3.0 ) { value = x * ( x + 3.0 ) * exp ( x ); } else { value = 2.0; } return value; } /******************************************************************************/ void test10 ( ) /******************************************************************************/ /* Purpose: TEST10 tests FEM1D_BVP_LINEAR. Discussion: We want to compute errors and do convergence rates for the following problem: - uxx + u = x for 0 < x < 1 u(0) = u(1) = 0 exact = x - sinh(x) / sinh(1) exact' = 1 - cosh(x) / sinh(1) Licensing: This code is distributed under the MIT license. Modified: 10 July 2015 Author: John Burkardt Reference: Dianne O'Leary, Scientific Computing with Case Studies, SIAM, 2008, ISBN13: 978-0-898716-66-5, LC: QA401.O44. */ { char command_filename[255]; FILE *command_unit; char data_filename[255]; FILE *data_unit; int e_log; int e_log_max = 6; double *h_plot; double h1; double *h1_plot; double l2; double *l2_plot; double mx; double *mx_plot; int n; int ne; int ne1; int ne2; int *ne_plot; char output_filename[255]; double r; double *u; double *x; double x_hi; double x_lo; printf ( "\n" ); printf ( "TEST10\n" ); printf ( " Solve -( A(x) U'(x) )' + C(x) U(x) = F(x)\n" ); printf ( " for 0 < x < 1, with U(0) = U(1) = 0.\n" ); printf ( " A(X) = 1.0\n" ); printf ( " C(X) = 1.0\n" ); printf ( " F(X) = X\n" ); printf ( " U(X) = X - SINH(X) / SINH(1)\n" ); printf ( "\n" ); printf ( " log(E) E L2error H1error Maxerror\n" ); printf ( "\n" ); h_plot = ( double * ) malloc ( ( e_log_max + 1 ) * sizeof ( double ) ); h1_plot = ( double * ) malloc ( ( e_log_max + 1 ) * sizeof ( double ) ); l2_plot = ( double * ) malloc ( ( e_log_max + 1 ) * sizeof ( double ) ); mx_plot = ( double * ) malloc ( ( e_log_max + 1 ) * sizeof ( double ) ); ne_plot = ( int * ) malloc ( ( e_log_max + 1 ) * sizeof ( int ) ); for ( e_log = 0; e_log <= e_log_max; e_log++ ) { ne = i4_power ( 2, e_log ); n = ne + 1; x_lo = 0.0; x_hi = 1.0; x = r8vec_linspace_new ( n, x_lo, x_hi ); u = fem1d_bvp_linear ( n, a10, c10, f10, x ); ne_plot[e_log] = ne; h_plot[e_log] = ( x_hi - x_lo ) / ( double ) ( ne ); l2 = l2_error_linear ( n, x, u, exact10 ); l2_plot[e_log] = l2; h1 = h1s_error_linear ( n, x, u, exact_ux10 ); h1_plot[e_log] = h1; mx = max_error_linear ( n, x, u, exact10 ); mx_plot[e_log] = mx; printf ( " %4d %4d %14.6g %14.6g %14.6g\n", e_log, ne, l2, h1, mx ); free ( u ); free ( x ); } printf ( "\n" ); printf ( " log(E1) E1 / E2 L2rate H1rate Maxrate\n" ); printf ( "\n" ); for ( e_log = 0; e_log < e_log_max; e_log++ ) { ne1 = ne_plot[e_log]; ne2 = ne_plot[e_log+1]; r = ( double ) ( ne2 ) / ( double ) ( ne1 ); l2 = l2_plot[e_log] / l2_plot[e_log+1]; l2 = log ( l2 ) / log ( r ); h1 = h1_plot[e_log] / h1_plot[e_log+1]; h1 = log ( h1 ) / log ( r ); mx = mx_plot[e_log] / mx_plot[e_log+1]; mx = log ( mx ) / log ( r ); printf ( " %4d %4d / %4d %14.6g %14.6g %14.6g\n", e_log, ne1, ne2, l2, h1, mx ); } /* Create the data file. */ strcpy ( data_filename, "data.txt" ); data_unit = fopen ( data_filename, "wt" ); for ( e_log = 0; e_log <= e_log_max; e_log++ ) { fprintf ( data_unit, " %d %g %g %g\n", ne_plot[e_log], l2_plot[e_log], h1_plot[e_log], mx_plot[e_log] ); } fclose ( data_unit ); printf ( "\n" ); printf ( " Created graphics data file \"%s\".\n", data_filename ); /* Plot the L2 error as a function of NE. */ strcpy ( command_filename, "commands_l2.txt" ); command_unit = fopen ( command_filename, "wt" ); strcpy ( output_filename, "l2.png" ); 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 '%s'\n", output_filename ); fprintf ( command_unit, "set xlabel '<---NE--->'\n" ); fprintf ( command_unit, "set ylabel '<---L2(NE)--->'\n" ); fprintf ( command_unit, "set title 'L2 error versus number of elements NE'\n" ); fprintf ( command_unit, "set logscale xy\n" ); fprintf ( command_unit, "set size ratio -1\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "plot '%s' using 1:2 with points pt 7 ps 2 lc rgb 'blue',\\\n", data_filename ); fprintf ( command_unit, " '%s' using 1:2 lw 3 linecolor rgb 'red'\n", data_filename ); fclose ( command_unit ); printf ( " Created graphics command file \"%s\".\n", command_filename ); /* Plot the H1 error as a function of NE. */ strcpy ( command_filename, "commands_h1.txt" ); command_unit = fopen ( command_filename, "wt" ); strcpy ( output_filename, "h1.png" ); 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 '%s'\n", output_filename ); fprintf ( command_unit, "set xlabel '<---NE--->'\n" ); fprintf ( command_unit, "set ylabel '<---H1(NE)--->'\n" ); fprintf ( command_unit, "set title 'H1 error versus number of elements NE'\n" ); fprintf ( command_unit, "set logscale xy\n" ); fprintf ( command_unit, "set size ratio -1\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "plot '%s' using 1:3 with points pt 7 ps 2 lc rgb 'blue',\\\n", data_filename ); fprintf ( command_unit, " '%s' using 1:3 lw 3 linecolor rgb 'red'\n", data_filename ); fclose ( command_unit ); printf ( " Created graphics command file \"%s\".\n", command_filename ); /* Plot the MX error as a function of NE. */ strcpy ( command_filename, "commands_mx.txt" ); command_unit = fopen ( command_filename, "wt" ); strcpy ( output_filename, "mx.png" ); 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 '%s'\n", output_filename ); fprintf ( command_unit, "set xlabel '<---NE--->'\n" ); fprintf ( command_unit, "set ylabel '<---MX(NE)--->'\n" ); fprintf ( command_unit, "set title 'Max error versus number of elements NE'\n" ); fprintf ( command_unit, "set logscale xy\n" ); fprintf ( command_unit, "set size ratio -1\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "plot '%s' using 1:4 with points pt 7 ps 2 lc rgb 'blue',\\\n", data_filename ); fprintf ( command_unit, " '%s' using 1:4 lw 3 linecolor rgb 'red'\n", data_filename ); fclose ( command_unit ); printf ( " Created graphics command file \"%s\".\n", command_filename ); /* Free memory. */ free ( h_plot ); free ( h1_plot ); free ( l2_plot ); free ( mx_plot ); free ( ne_plot ); return; } /******************************************************************************/ double a10 ( double x ) /******************************************************************************/ /* Purpose: A10 evaluates A function #10. Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double A10, the value of A(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double c10 ( double x ) /******************************************************************************/ /* Purpose: C10 evaluates C function #10. Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double C10, the value of C(X). */ { double value; value = 1.0; return value; } /******************************************************************************/ double exact10 ( double x ) /******************************************************************************/ /* Purpose: EXACT10 evaluates exact solution #10. Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT10, the value of U(X). */ { double value; value = x - sinh ( x ) / sinh ( 1.0 ); return value; } /******************************************************************************/ double exact_ux10 ( double x ) /******************************************************************************/ /* Purpose: EXACT_UX10 evaluates the derivative of exact solution #10. Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double EXACT_UX10, the value of dUdX(X). */ { double value; value = 1.0 - cosh ( x ) / sinh ( 1.0 ); return value; } /******************************************************************************/ double f10 ( double x ) /******************************************************************************/ /* Purpose: F10 evaluates right hand side function #10. Licensing: This code is distributed under the MIT license. Modified: 18 July 2015 Author: John Burkardt Parameters: Input, double X, the evaluation point. Output, double F10, the value of F(X). */ { double value; value = x; return value; }