Last active
November 12, 2021 13:12
-
-
Save mdecourse/b5df54ab3d5f2f079f785541d1178a66 to your computer and use it in GitHub Desktop.
Optimization Programs
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # include <math.h> | |
| # include <stdlib.h> | |
| # include <stdio.h> | |
| # include <time.h> | |
| //# include "fem1d_bvp_linear.h" | |
| double *fem1d_bvp_linear ( int n, double a ( double x ), double c ( double x ), | |
| double f ( double x ), double x[] ); | |
| double h1s_error_linear ( int n, double x[], double u[], | |
| double exact_ux ( double x ) ); | |
| int i4_power ( int i, int j ); | |
| int *i4vec_zero_new ( int n ); | |
| double l1_error ( int n, double x[], double u[], | |
| double exact ( double x ) ); | |
| double l2_error_linear ( int n, double x[], double u[], | |
| double exact ( double x ) ); | |
| double max_error_linear ( int n, double x[], double u[], | |
| double exact ( double x ) ); | |
| double r8_max ( double x, double y ); | |
| double *r8mat_solve2 ( int n, double a[], double b[], int *ierror ); | |
| double *r8mat_zero_new ( int m, int n ); | |
| double *r8vec_linspace_new ( int n, double alo, double ahi ); | |
| double *r8vec_zero_new ( int n ); | |
| void timestamp ( ); | |
| /******************************************************************************/ | |
| double *fem1d_bvp_linear ( int n, double a ( double x ), double c ( double x ), | |
| double f ( double x ), double x[] ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| FEM1D_BVP_LINEAR solves a two point boundary value problem. | |
| Location: | |
| http://people.sc.fsu.edu/~jburkardt/c_src/fem1d_bvp_linear/fem1d_bvp_linear.c | |
| Discussion: | |
| The program uses the finite element method, with piecewise linear basis | |
| functions to solve a boundary value problem in one dimension. | |
| The problem is defined on the region 0 <= x <= 1. | |
| The following differential equation is imposed between 0 and 1: | |
| - d/dx a(x) du/dx + c(x) * u(x) = f(x) | |
| where a(x), c(x), and f(x) are given functions. | |
| At the boundaries, the following conditions are applied: | |
| u(0.0) = 0.0 | |
| u(1.0) = 0.0 | |
| A set of N equally spaced nodes is defined on this | |
| interval, with 0 = X(1) < X(2) < ... < X(N) = 1.0. | |
| At each node I, we associate a piecewise linear basis function V(I,X), | |
| which is 0 at all nodes except node I. This implies that V(I,X) is | |
| everywhere 0 except that | |
| for X(I-1) <= X <= X(I): | |
| V(I,X) = ( X - X(I-1) ) / ( X(I) - X(I-1) ) | |
| for X(I) <= X <= X(I+1): | |
| V(I,X) = ( X(I+1) - X ) / ( X(I+1) - X(I) ) | |
| We now assume that the solution U(X) can be written as a linear | |
| sum of these basis functions: | |
| U(X) = sum ( 1 <= J <= N ) U(J) * V(J,X) | |
| where U(X) on the left is the function of X, but on the right, | |
| is meant to indicate the coefficients of the basis functions. | |
| To determine the coefficient U(J), we multiply the original | |
| differential equation by the basis function V(J,X), and use | |
| integration by parts, to arrive at the I-th finite element equation: | |
| Integral A(X) * U'(X) * V'(I,X) + C(X) * U(X) * V(I,X) dx | |
| = Integral F(X) * V(I,X) dx | |
| We note that the functions U(X) and U'(X) can be replaced by | |
| the finite element form involving the linear sum of basis functions, | |
| but we also note that the resulting integrand will only be nonzero | |
| for terms where J = I - 1, I, or I + 1. | |
| By writing this equation for basis functions I = 2 through N - 1, | |
| and using the boundary conditions, we have N linear equations | |
| for the N unknown coefficients U(1) through U(N), which can | |
| be easily solved. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 18 June 2014 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of nodes. | |
| Input, double A ( double X ), evaluates a(x); | |
| Input, double C ( double X ), evaluates c(x); | |
| Input, double F ( double X ), evaluates f(x); | |
| Input, double X[N], the mesh points. | |
| Output, double FEM1D_BVP_LINEAR[N], the finite element coefficients, | |
| which are also the value of the computed solution at the mesh points. | |
| */ | |
| { | |
| # define QUAD_NUM 2 | |
| double abscissa[QUAD_NUM] = { | |
| -0.577350269189625764509148780502, | |
| +0.577350269189625764509148780502 }; | |
| double *amat; | |
| double axq; | |
| double *b; | |
| double cxq; | |
| int e; | |
| int e_num; | |
| double fxq; | |
| int i; | |
| int ierror; | |
| int j; | |
| int l; | |
| int q; | |
| int quad_num = QUAD_NUM; | |
| int r; | |
| double *u; | |
| double weight[QUAD_NUM] = { 1.0, 1.0 }; | |
| double wq; | |
| double vl; | |
| double vlp; | |
| double vr; | |
| double vrp; | |
| double xl; | |
| double xq; | |
| double xr; | |
| /* | |
| Zero out the matrix and right hand side. | |
| */ | |
| amat = r8mat_zero_new ( n, n ); | |
| b = r8vec_zero_new ( n ); | |
| e_num = n - 1; | |
| for ( e = 0; e < e_num; e++ ) | |
| { | |
| l = e; | |
| r = e + 1; | |
| xl = x[l]; | |
| xr = x[r]; | |
| for ( q = 0; q < quad_num; q++ ) | |
| { | |
| xq = ( ( 1.0 - abscissa[q] ) * xl | |
| + ( 1.0 + abscissa[q] ) * xr ) | |
| / 2.0; | |
| wq = weight[q] * ( xr - xl ) / 2.0; | |
| vl = ( xr - xq ) / ( xr - xl ); | |
| vlp = - 1.0 / ( xr - xl ); | |
| vr = ( xq - xl ) / ( xr - xl ); | |
| vrp = + 1.0 / ( xr - xl ); | |
| axq = a ( xq ); | |
| cxq = c ( xq ); | |
| fxq = f ( xq ); | |
| amat[l+l*n] = amat[l+l*n] + wq * ( vlp * axq * vlp + vl * cxq * vl ); | |
| amat[l+r*n] = amat[l+r*n] + wq * ( vlp * axq * vrp + vl * cxq * vr ); | |
| b[l] = b[l] + wq * ( vl * fxq ); | |
| amat[r+l*n] = amat[r+l*n] + wq * ( vrp * axq * vlp + vr * cxq * vl ); | |
| amat[r+r*n] = amat[r+r*n] + wq * ( vrp * axq * vrp + vr * cxq * vr ); | |
| b[r] = b[r] + wq * ( vr * fxq ); | |
| } | |
| } | |
| /* | |
| Equation 1 is the left boundary condition, U(0.0) = 0.0; | |
| */ | |
| for ( j = 0; j < n; j++ ) | |
| { | |
| amat[0+j*n] = 0.0; | |
| } | |
| b[0] = 0.0; | |
| for ( i = 1; i < n; i++ ) | |
| { | |
| b[i] = b[i] - amat[i+0*n] * b[0]; | |
| } | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| amat[i+0*n] = 0.0; | |
| } | |
| amat[0+0*n] = 1.0; | |
| /* | |
| Equation N is the right boundary condition, U(1.0) = 0.0; | |
| */ | |
| for ( j = 0; j < n; j++ ) | |
| { | |
| amat[n-1+j*n] = 0.0; | |
| } | |
| b[n-1] = 0.0; | |
| for ( i = 0; i < n - 1; i++ ) | |
| { | |
| b[i] = b[i] - amat[i+(n-1)*n] * b[n-1]; | |
| } | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| amat[i+(n-1)*n] = 0.0; | |
| } | |
| amat[n-1+(n-1)*n] = 1.0; | |
| /* | |
| Solve the linear system. | |
| */ | |
| u = r8mat_solve2 ( n, amat, b, &ierror ); | |
| free ( amat ); | |
| free ( b ); | |
| return u; | |
| # undef QUAD_NUM | |
| } | |
| /******************************************************************************/ | |
| double h1s_error_linear ( int n, double x[], double u[], | |
| double exact_ux ( double x ) ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| H1S_ERROR_LINEAR estimates the seminorm error of a finite element solution. | |
| Discussion: | |
| We assume the finite element method has been used, over an interval [A,B] | |
| involving N nodes, with piecewise linear elements used for the basis. | |
| The coefficients U(1:N) have been computed, and a formula for the | |
| exact derivative is known. | |
| This function estimates the seminorm of the error: | |
| SEMINORM = Integral ( A <= X <= B ) ( dU(X)/dx - EXACT_UX(X) )^2 dX | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 18 June 2014 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of nodes. | |
| Input, double X(N), the mesh points. | |
| Input, double U(N), the finite element coefficients. | |
| Input, function EQ = EXACT_UX ( X ), returns the value of the exact | |
| derivative at the point X. | |
| Output, double H1S_ERROR_LINEAR, the estimated seminorm of | |
| the error. | |
| */ | |
| { | |
| # define QUAD_NUM 2 | |
| double abscissa[QUAD_NUM] = { | |
| -0.577350269189625764509148780502, | |
| +0.577350269189625764509148780502 }; | |
| double exq; | |
| double h1s; | |
| int i; | |
| int q; | |
| int quad_num = QUAD_NUM; | |
| double ul; | |
| double ur; | |
| double uxq; | |
| double weight[QUAD_NUM] = { 1.0, 1.0 }; | |
| double wq; | |
| double xl; | |
| double xq; | |
| double xr; | |
| h1s = 0.0; | |
| /* | |
| Integrate over each interval. | |
| */ | |
| for ( i = 0; i < n - 1; i++ ) | |
| { | |
| xl = x[i]; | |
| xr = x[i+1]; | |
| ul = u[i]; | |
| ur = u[i+1]; | |
| for ( q = 0; q < quad_num; q++ ) | |
| { | |
| xq = ( ( 1.0 - abscissa[q] ) * xl | |
| + ( 1.0 + abscissa[q] ) * xr ) | |
| / 2.0; | |
| wq = weight[q] * ( xr - xl ) / 2.0; | |
| /* | |
| The piecewise linear derivative is a constant in the interval. | |
| */ | |
| uxq = ( ur - ul ) / ( xr - xl ); | |
| exq = exact_ux ( xq ); | |
| h1s = h1s + wq * pow ( uxq - exq, 2); | |
| } | |
| } | |
| h1s = sqrt ( h1s ); | |
| return h1s; | |
| # undef QUAD_NUM | |
| } | |
| /******************************************************************************/ | |
| int i4_power ( int i, int j ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| I4_POWER returns the value of I^J. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 23 October 2007 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int I, J, the base and the power. J should be nonnegative. | |
| Output, int I4_POWER, the value of I^J. | |
| */ | |
| { | |
| int k; | |
| int value; | |
| if ( j < 0 ) | |
| { | |
| if ( i == 1 ) | |
| { | |
| value = 1; | |
| } | |
| else if ( i == 0 ) | |
| { | |
| fprintf ( stderr, "\n" ); | |
| fprintf ( stderr, "I4_POWER - Fatal error!\n" ); | |
| fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" ); | |
| exit ( 1 ); | |
| } | |
| else | |
| { | |
| value = 0; | |
| } | |
| } | |
| else if ( j == 0 ) | |
| { | |
| if ( i == 0 ) | |
| { | |
| fprintf ( stderr, "\n" ); | |
| fprintf ( stderr, "I4_POWER - Fatal error!\n" ); | |
| fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" ); | |
| exit ( 1 ); | |
| } | |
| else | |
| { | |
| value = 1; | |
| } | |
| } | |
| else if ( j == 1 ) | |
| { | |
| value = i; | |
| } | |
| else | |
| { | |
| value = 1; | |
| for ( k = 1; k <= j; k++ ) | |
| { | |
| value = value * i; | |
| } | |
| } | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| int *i4vec_zero_new ( int n ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| I4VEC_ZERO_NEW creates and zeroes an I4VEC. | |
| Discussion: | |
| An I4VEC is a vector of I4's. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 05 September 2008 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of entries in the vector. | |
| Output, int I4VEC_ZERO_NEW[N], a vector of zeroes. | |
| */ | |
| { | |
| int *a; | |
| int i; | |
| a = ( int * ) malloc ( n * sizeof ( int ) ); | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| a[i] = 0; | |
| } | |
| return a; | |
| } | |
| /******************************************************************************/ | |
| double l1_error ( int n, double x[], double u[], double exact ( double x ) ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| L1_ERROR estimates the l1 error norm of a finite element solution. | |
| Discussion: | |
| We assume the finite element method has been used, over an interval [A,B] | |
| involving N nodes. | |
| The coefficients U(1:N) have been computed, and a formula for the | |
| exact solution is known. | |
| This function estimates the little l1 norm of the error: | |
| L1_NORM = sum ( 1 <= I <= N ) abs ( U(i) - EXACT(X(i)) ) | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 14 June 2014 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of nodes. | |
| Input, double X[N], the mesh points. | |
| Input, double U[N], the finite element coefficients. | |
| Input, function EQ = EXACT ( X ), returns the value of the exact | |
| solution at the point X. | |
| Output, double L1_ERROR, the little l1 norm of the error. | |
| */ | |
| { | |
| int i; | |
| double e1; | |
| e1 = 0.0; | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| e1 = e1 + fabs ( u[i] - exact ( x[i] ) ); | |
| } | |
| e1 = e1 / ( double ) n; | |
| return e1; | |
| } | |
| /******************************************************************************/ | |
| double l2_error_linear ( int n, double x[], double u[], | |
| double exact ( double x ) ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| L2_ERROR_LINEAR estimates the L2 error norm of a finite element solution. | |
| Discussion: | |
| We assume the finite element method has been used, over an interval [A,B] | |
| involving N nodes, with piecewise linear elements used for the basis. | |
| The coefficients U(1:N) have been computed, and a formula for the | |
| exact solution is known. | |
| This function estimates the L2 norm of the error: | |
| L2_NORM = Integral ( A <= X <= B ) ( U(X) - EXACT(X) )^2 dX | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 18 June 2014 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of nodes. | |
| Input, double X[N], the mesh points. | |
| Input, double U[N], the finite element coefficients. | |
| Input, function EQ = EXACT ( X ), returns the value of the exact | |
| solution at the point X. | |
| Output, double L2_ERROR_LINEAR, the estimated L2 norm of the error. | |
| */ | |
| { | |
| # define QUAD_NUM 2 | |
| double abscissa[QUAD_NUM] = { | |
| -0.577350269189625764509148780502, | |
| +0.577350269189625764509148780502 }; | |
| double e2; | |
| double eq; | |
| int i; | |
| int q; | |
| int quad_num = QUAD_NUM; | |
| double ul; | |
| double ur; | |
| double uq; | |
| double weight[QUAD_NUM] = { 1.0, 1.0 }; | |
| double wq; | |
| double xl; | |
| double xq; | |
| double xr; | |
| e2 = 0.0; | |
| /* | |
| Integrate over each interval. | |
| */ | |
| for ( i = 0; i < n - 1; i++ ) | |
| { | |
| xl = x[i]; | |
| xr = x[i+1]; | |
| ul = u[i]; | |
| ur = u[i+1]; | |
| for ( q = 0; q < quad_num; q++ ) | |
| { | |
| xq = ( ( 1.0 - abscissa[q] ) * xl | |
| + ( 1.0 + abscissa[q] ) * xr ) | |
| / 2.0; | |
| wq = weight[q] * ( xr - xl ) / 2.0; | |
| /* | |
| Use the fact that U is a linear combination of piecewise linears. | |
| */ | |
| uq = ( ( xr - xq ) * ul | |
| + ( xq - xl ) * ur ) | |
| / ( xr - xl ); | |
| eq = exact ( xq ); | |
| e2 = e2 + wq * pow ( uq - eq, 2 ); | |
| } | |
| } | |
| e2 = sqrt ( e2 ); | |
| return e2; | |
| # undef QUAD_NUM | |
| } | |
| /******************************************************************************/ | |
| double max_error_linear ( int n, double x[], double u[], | |
| double exact ( double x ) ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| MAX_ERROR_LINEAR estimates the max error norm of a finite element solution. | |
| Discussion: | |
| We assume the finite element method has been used, over an interval [A,B] | |
| involving N nodes, with piecewise linear elements used for the basis. | |
| The coefficients U(1:N) have been computed, and a formula for the | |
| exact solution is known. | |
| This function estimates the max norm of the error: | |
| MAX_NORM = Integral ( A <= X <= B ) max ( abs ( U(X) - EXACT(X) ) ) dX | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 08 July 2015 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of nodes. | |
| Input, double X[N], the mesh points. | |
| Input, double U[N], the finite element coefficients. | |
| Input, function EQ = EXACT ( X ), returns the value of the exact | |
| solution at the point X. | |
| Output, double MAX_ERROR_LINEAR, the estimated max norm of the error. | |
| */ | |
| { | |
| int e; | |
| int e_num; | |
| double eq; | |
| int l; | |
| int q; | |
| int quad_num = 8; | |
| int r; | |
| double ul; | |
| double ur; | |
| double uq; | |
| double value; | |
| double xl; | |
| double xq; | |
| double xr; | |
| value = 0.0; | |
| /* | |
| Integrate over each interval. | |
| */ | |
| e_num = n - 1; | |
| for ( e = 0; e < e_num; e++ ) | |
| { | |
| l = e; | |
| xl = x[l]; | |
| ul = u[l]; | |
| r = e + 1; | |
| xr = x[r]; | |
| ur = u[r]; | |
| for ( q = 0; q < quad_num; q++ ) | |
| { | |
| xq = ( ( double ) ( quad_num - q ) * xl | |
| + ( double ) ( q ) * xr ) | |
| / ( double ) ( quad_num ); | |
| /* | |
| Use the fact that U is a linear combination of piecewise linears. | |
| */ | |
| uq = ( ( xr - xq ) * ul | |
| + ( xq - xl ) * ur ) | |
| / ( xr - xl ); | |
| eq = exact ( xq ); | |
| value = r8_max ( value, fabs ( uq - eq ) ); | |
| } | |
| } | |
| /* | |
| For completeness, check last node. | |
| */ | |
| xq = x[n-1]; | |
| uq = u[n-1]; | |
| eq = exact ( xq ); | |
| value = r8_max ( value, fabs ( uq - eq ) ); | |
| /* | |
| Integral approximation requires multiplication by interval length. | |
| */ | |
| value = value * ( x[n-1] - x[0] ); | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| double r8_max ( double x, double y ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| R8_MAX returns the maximum of two R8's. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 07 May 2006 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, double X, Y, the quantities to compare. | |
| Output, double R8_MAX, the maximum of X and Y. | |
| */ | |
| { | |
| double value; | |
| if ( y < x ) | |
| { | |
| value = x; | |
| } | |
| else | |
| { | |
| value = y; | |
| } | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| double *r8mat_solve2 ( int n, double a[], double b[], int *ierror ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| R8MAT_SOLVE2 computes the solution of an N by N linear system. | |
| Discussion: | |
| An R8MAT is a doubly dimensioned array of R8 values, stored as a vector | |
| in column-major order. | |
| The linear system may be represented as | |
| A*X = B | |
| If the linear system is singular, but consistent, then the routine will | |
| still produce a solution. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 20 August 2010 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of equations. | |
| Input/output, double A[N*N]. | |
| On input, A is the coefficient matrix to be inverted. | |
| On output, A has been overwritten. | |
| Input/output, double B[N]. | |
| On input, B is the right hand side of the system. | |
| On output, B has been overwritten. | |
| Output, double R8MAT_SOLVE2[N], the solution of the linear system. | |
| Output, int *IERROR. | |
| 0, no error detected. | |
| 1, consistent singularity. | |
| 2, inconsistent singularity. | |
| */ | |
| { | |
| double amax; | |
| int i; | |
| int imax; | |
| int j; | |
| int k; | |
| int *piv; | |
| double *x; | |
| *ierror = 0; | |
| piv = i4vec_zero_new ( n ); | |
| x = r8vec_zero_new ( n ); | |
| /* | |
| Process the matrix. | |
| */ | |
| for ( k = 1; k <= n; k++ ) | |
| { | |
| /* | |
| In column K: | |
| Seek the row IMAX with the properties that: | |
| IMAX has not already been used as a pivot; | |
| A(IMAX,K) is larger in magnitude than any other candidate. | |
| */ | |
| amax = 0.0; | |
| imax = 0; | |
| for ( i = 1; i <= n; i++ ) | |
| { | |
| if ( piv[i-1] == 0 ) | |
| { | |
| if ( amax < fabs ( a[i-1+(k-1)*n] ) ) | |
| { | |
| imax = i; | |
| amax = fabs ( a[i-1+(k-1)*n] ); | |
| } | |
| } | |
| } | |
| /* | |
| If you found a pivot row IMAX, then, | |
| eliminate the K-th entry in all rows that have not been used for pivoting. | |
| */ | |
| if ( imax != 0 ) | |
| { | |
| piv[imax-1] = k; | |
| for ( j = k+1; j <= n; j++ ) | |
| { | |
| a[imax-1+(j-1)*n] = a[imax-1+(j-1)*n] / a[imax-1+(k-1)*n]; | |
| } | |
| b[imax-1] = b[imax-1] / a[imax-1+(k-1)*n]; | |
| a[imax-1+(k-1)*n] = 1.0; | |
| for ( i = 1; i <= n; i++ ) | |
| { | |
| if ( piv[i-1] == 0 ) | |
| { | |
| for ( j = k+1; j <= n; j++ ) | |
| { | |
| a[i-1+(j-1)*n] = a[i-1+(j-1)*n] - a[i-1+(k-1)*n] * a[imax-1+(j-1)*n]; | |
| } | |
| b[i-1] = b[i-1] - a[i-1+(k-1)*n] * b[imax-1]; | |
| a[i-1+(k-1)*n] = 0.0; | |
| } | |
| } | |
| } | |
| } | |
| /* | |
| Now, every row with nonzero PIV begins with a 1, and | |
| all other rows are all zero. Begin solution. | |
| */ | |
| for ( j = n; 1 <= j; j-- ) | |
| { | |
| imax = 0; | |
| for ( k = 1; k <= n; k++ ) | |
| { | |
| if ( piv[k-1] == j ) | |
| { | |
| imax = k; | |
| } | |
| } | |
| if ( imax == 0 ) | |
| { | |
| x[j-1] = 0.0; | |
| if ( b[j-1] == 0.0 ) | |
| { | |
| *ierror = 1; | |
| printf ( "\n" ); | |
| printf ( "R8MAT_SOLVE2 - Warning:\n" ); | |
| printf ( " Consistent singularity, equation = %d\n", j ); | |
| } | |
| else | |
| { | |
| *ierror = 2; | |
| printf ( "\n" ); | |
| printf ( "R8MAT_SOLVE2 - Warning:\n" ); | |
| printf ( " Inconsistent singularity, equation = %d\n", j ); | |
| } | |
| } | |
| else | |
| { | |
| x[j-1] = b[imax-1]; | |
| for ( i = 1; i <= n; i++ ) | |
| { | |
| if ( i != imax ) | |
| { | |
| b[i-1] = b[i-1] - a[i-1+(j-1)*n] * x[j-1]; | |
| } | |
| } | |
| } | |
| } | |
| free ( piv ); | |
| return x; | |
| } | |
| /******************************************************************************/ | |
| double *r8mat_zero_new ( int m, int n ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| R8MAT_ZERO_NEW returns a new zeroed R8MAT. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 26 September 2008 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int M, N, the number of rows and columns. | |
| Output, double R8MAT_ZERO[M*N], the new zeroed matrix. | |
| */ | |
| { | |
| double *a; | |
| int i; | |
| int j; | |
| a = ( double * ) malloc ( m * n * sizeof ( double ) ); | |
| for ( j = 0; j < n; j++ ) | |
| { | |
| for ( i = 0; i < m; i++ ) | |
| { | |
| a[i+j*m] = 0.0; | |
| } | |
| } | |
| return a; | |
| } | |
| /******************************************************************************/ | |
| double *r8vec_linspace_new ( int n, double a, double b ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| R8VEC_LINSPACE_NEW creates a vector of linearly spaced values. | |
| Discussion: | |
| An R8VEC is a vector of R8's. | |
| 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. | |
| In other words, the interval is divided into N-1 even subintervals, | |
| and the endpoints of intervals are used as the points. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 March 2011 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of entries in the vector. | |
| Input, double A, B, the first and last entries. | |
| Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. | |
| */ | |
| { | |
| int i; | |
| double *x; | |
| x = ( double * ) malloc ( n * sizeof ( double ) ); | |
| if ( n == 1 ) | |
| { | |
| x[0] = ( a + b ) / 2.0; | |
| } | |
| else | |
| { | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| x[i] = ( ( double ) ( n - 1 - i ) * a | |
| + ( double ) ( i ) * b ) | |
| / ( double ) ( n - 1 ); | |
| } | |
| } | |
| return x; | |
| } | |
| /******************************************************************************/ | |
| double *r8vec_zero_new ( int n ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| R8VEC_ZERO_NEW creates and zeroes an R8VEC. | |
| Discussion: | |
| An R8VEC is a vector of R8's. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 25 March 2009 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, int N, the number of entries in the vector. | |
| Output, double R8VEC_ZERO_NEW[N], a vector of zeroes. | |
| */ | |
| { | |
| double *a; | |
| int i; | |
| a = ( double * ) malloc ( n * sizeof ( double ) ); | |
| for ( i = 0; i < n; i++ ) | |
| { | |
| a[i] = 0.0; | |
| } | |
| return a; | |
| } | |
| /******************************************************************************/ | |
| void timestamp ( void ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| TIMESTAMP prints the current YMDHMS date as a time stamp. | |
| Example: | |
| 31 May 2001 09:45:54 AM | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 24 September 2003 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| None | |
| */ | |
| { | |
| # 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 ); | |
| fprintf ( stdout, "%s\n", time_buffer ); | |
| return; | |
| # undef TIME_SIZE | |
| } | |
| //####################### | |
| //####################### | |
| # include <stdlib.h> | |
| # include <stdio.h> | |
| # include <time.h> | |
| int main ( void ); | |
| void assemble ( double adiag[], double aleft[], double arite[], double f[], | |
| double h[], int indx[], int nl, int node[], int nu, int nquad, int nsub, | |
| double ul, double ur, double xn[], double xquad[] ); | |
| double ff ( double x ); | |
| void geometry ( double h[], int ibc, int indx[], int nl, int node[], int nsub, | |
| int *nu, double xl, double xn[], double xquad[], double xr ); | |
| void init ( int *ibc, int *nquad, double *ul, double *ur, double *xl, | |
| double *xr ); | |
| void output ( double f[], int ibc, int indx[], int nsub, int nu, double ul, | |
| double ur, double xn[] ); | |
| void phi ( int il, double x, double *phii, double *phiix, double xleft, | |
| double xrite ); | |
| double pp ( double x ); | |
| void prsys ( double adiag[], double aleft[], double arite[], double f[], | |
| int nu ); | |
| double qq ( double x ); | |
| void solve ( double adiag[], double aleft[], double arite[], double f[], | |
| int nu ); | |
| void timestamp ( ); | |
| /******************************************************************************/ | |
| int main ( ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| MAIN is the main program for FEM1D. | |
| Discussion: | |
| FEM1D solves a one dimensional ODE using the finite element method. | |
| The differential equation solved is | |
| - d/dX (P dU/dX) + Q U = F | |
| The finite-element method uses piecewise linear basis functions. | |
| Here U is an unknown scalar function of X defined on the | |
| interval [XL,XR], and P, Q and F are given functions of X. | |
| The values of U or U' at XL and XR are also specified. | |
| The interval [XL,XR] is "meshed" with NSUB+1 points, | |
| XN(0) = XL, XN(1)=XL+H, XN(2)=XL+2*H, ..., XN(NSUB)=XR. | |
| This creates NSUB subintervals, with interval number 1 | |
| having endpoints XN(0) and XN(1), and so on up to interval | |
| NSUB, which has endpoints XN(NSUB-1) and XN(NSUB). | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| double ADIAG(NU), the "diagonal" coefficients. That is, ADIAG(I) is | |
| the coefficient of the I-th unknown in the I-th equation. | |
| double ALEFT(NU), the "left hand" coefficients. That is, ALEFT(I) | |
| is the coefficient of the (I-1)-th unknown in the I-th equation. | |
| There is no value in ALEFT(1), since the first equation | |
| does not refer to a "0-th" unknown. | |
| double ARITE(NU). | |
| ARITE(I) is the "right hand" coefficient of the I-th | |
| equation in the linear system. ARITE(I) is the coefficient | |
| of the (I+1)-th unknown in the I-th equation. There is | |
| no value in ARITE(NU) because the NU-th equation does not | |
| refer to an "NU+1"-th unknown. | |
| double F(NU). | |
| ASSEMBLE stores into F the right hand side of the linear | |
| equations. | |
| SOLVE replaces those values of F by the solution of the | |
| linear equations. | |
| double H(NSUB) | |
| H(I) is the length of subinterval I. This code uses | |
| equal spacing for all the subintervals. | |
| int IBC. | |
| IBC declares what the boundary conditions are. | |
| 1, at the left endpoint, U has the value UL, | |
| at the right endpoint, U' has the value UR. | |
| 2, at the left endpoint, U' has the value UL, | |
| at the right endpoint, U has the value UR. | |
| 3, at the left endpoint, U has the value UL, | |
| and at the right endpoint, U has the value UR. | |
| 4, at the left endpoint, U' has the value UL, | |
| at the right endpoint U' has the value UR. | |
| int INDX[NSUB+1]. | |
| For a node I, INDX(I) is the index of the unknown | |
| associated with node I. | |
| If INDX(I) is equal to -1, then no unknown is associated | |
| with the node, because a boundary condition fixing the | |
| value of U has been applied at the node instead. | |
| Unknowns are numbered beginning with 1. | |
| If IBC is 2 or 4, then there is an unknown value of U | |
| at node 0, which will be unknown number 1. Otherwise, | |
| unknown number 1 will be associated with node 1. | |
| If IBC is 1 or 4, then there is an unknown value of U | |
| at node NSUB, which will be unknown NSUB or NSUB+1, | |
| depending on whether there was an unknown at node 0. | |
| int NL. | |
| The number of basis functions used in a single | |
| subinterval. (NL-1) is the degree of the polynomials | |
| used. For this code, NL is fixed at 2, meaning that | |
| piecewise linear functions are used as the basis. | |
| int NODE[NL*NSUB]. | |
| For each subinterval I: | |
| NODE[0+I*2] is the number of the left node, and | |
| NODE[1+I*2] is the number of the right node. | |
| int NQUAD. | |
| The number of quadrature points used in a subinterval. | |
| This code uses NQUAD = 1. | |
| int NSUB. | |
| The number of subintervals into which the interval [XL,XR] is broken. | |
| int NU. | |
| NU is the number of unknowns in the linear system. | |
| Depending on the value of IBC, there will be NSUB-1, | |
| NSUB, or NSUB+1 unknown values, which are the coefficients | |
| of basis functions. | |
| double UL. | |
| If IBC is 1 or 3, UL is the value that U is required | |
| to have at X = XL. | |
| If IBC is 2 or 4, UL is the value that U' is required | |
| to have at X = XL. | |
| double UR. | |
| If IBC is 2 or 3, UR is the value that U is required | |
| to have at X = XR. | |
| If IBC is 1 or 4, UR is the value that U' is required | |
| to have at X = XR. | |
| double XL. | |
| XL is the left endpoint of the interval over which the | |
| differential equation is being solved. | |
| double XN(0:NSUB). | |
| XN(I) is the location of the I-th node. XN(0) is XL, | |
| and XN(NSUB) is XR. | |
| double XQUAD(NSUB) | |
| XQUAD(I) is the location of the single quadrature point | |
| in interval I. | |
| double XR. | |
| XR is the right endpoint of the interval over which the | |
| differential equation is being solved. | |
| */ | |
| { | |
| # define NSUB 5 | |
| # define NL 2 | |
| double adiag[NSUB+1]; | |
| double aleft[NSUB+1]; | |
| double arite[NSUB+1]; | |
| double f[NSUB+1]; | |
| double h[NSUB]; | |
| int ibc; | |
| int indx[NSUB+1]; | |
| int node[NL*NSUB]; | |
| int nquad; | |
| int nu; | |
| double ul; | |
| double ur; | |
| double xl; | |
| double xn[NSUB+1]; | |
| double xquad[NSUB]; | |
| double xr; | |
| timestamp ( ); | |
| printf ( "\n" ); | |
| printf ( "FEM1D\n" ); | |
| printf ( " C version\n" ); | |
| printf ( "\n" ); | |
| printf ( " Solve the two-point boundary value problem\n" ); | |
| printf ( "\n" ); | |
| printf ( " - d/dX (P dU/dX) + Q U = F\n" ); | |
| printf ( "\n" ); | |
| printf ( " on the interval [XL,XR], specifying\n" ); | |
| printf ( " the value of U or U' at each end.\n" ); | |
| printf ( "\n" ); | |
| printf ( " The interval [XL,XR] is broken into NSUB = %d subintervals\n", NSUB ); | |
| printf ( " Number of basis functions per element is NL = %d\n", NL ); | |
| /* | |
| Initialize the data. | |
| */ | |
| init ( &ibc, &nquad, &ul, &ur, &xl, &xr ); | |
| /* | |
| Compute the geometric quantities. | |
| */ | |
| geometry ( h, ibc, indx, NL, node, NSUB, &nu, xl, xn, xquad, xr ); | |
| /* | |
| Assemble the linear system. | |
| */ | |
| assemble ( adiag, aleft, arite, f, h, indx, NL, node, nu, nquad, | |
| NSUB, ul, ur, xn, xquad ); | |
| /* | |
| Print out the linear system. | |
| */ | |
| prsys ( adiag, aleft, arite, f, nu ); | |
| /* | |
| Solve the linear system. | |
| */ | |
| solve ( adiag, aleft, arite, f, nu ); | |
| /* | |
| Print out the solution. | |
| */ | |
| output ( f, ibc, indx, NSUB, nu, ul, ur, xn ); | |
| /* | |
| Terminate. | |
| */ | |
| printf ( "\n" ); | |
| printf ( "FEM1D:\n" ); | |
| printf ( " Normal end of execution.\n" ); | |
| printf ( "\n" ); | |
| timestamp ( ); | |
| return 0; | |
| # undef NL | |
| # undef NSUB | |
| } | |
| /******************************************************************************/ | |
| void assemble ( double adiag[], double aleft[], double arite[], double f[], | |
| double h[], int indx[], int nl, int node[], int nu, int nquad, int nsub, | |
| double ul, double ur, double xn[], double xquad[] ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| ASSEMBLE assembles the matrix and right-hand-side of the linear system. | |
| Discussion: | |
| The linear system has the form: | |
| K * C = F | |
| that is to be solved for the coefficients C. | |
| Numerical integration is used to compute the entries of K and F. | |
| Note that a 1 point quadrature rule, which is sometimes used to | |
| assemble the matrix and right hand side, is just barely accurate | |
| enough for simple problems. If you want better results, you | |
| should use a quadrature rule that is more accurate. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Output, double ADIAG(NU), the "diagonal" coefficients. That is, | |
| ADIAG(I) is the coefficient of the I-th unknown in the I-th equation. | |
| Output, double ALEFT(NU), the "left hand" coefficients. That is, | |
| ALEFT(I) is the coefficient of the (I-1)-th unknown in the I-th equation. | |
| There is no value in ALEFT(1), since the first equation | |
| does not refer to a "0-th" unknown. | |
| Output, double ARITE(NU). | |
| ARITE(I) is the "right hand" coefficient of the I-th | |
| equation in the linear system. ARITE(I) is the coefficient | |
| of the (I+1)-th unknown in the I-th equation. There is | |
| no value in ARITE(NU) because the NU-th equation does not | |
| refer to an "NU+1"-th unknown. | |
| Output, double F(NU). | |
| ASSEMBLE stores into F the right hand side of the linear | |
| equations. | |
| SOLVE replaces those values of F by the solution of the | |
| linear equations. | |
| Input, double H(NSUB) | |
| H(I) is the length of subinterval I. This code uses | |
| equal spacing for all the subintervals. | |
| Input, int INDX[NSUB+1]. | |
| For a node I, INDX(I) is the index of the unknown | |
| associated with node I. | |
| If INDX(I) is equal to -1, then no unknown is associated | |
| with the node, because a boundary condition fixing the | |
| value of U has been applied at the node instead. | |
| Unknowns are numbered beginning with 1. | |
| If IBC is 2 or 4, then there is an unknown value of U | |
| at node 0, which will be unknown number 1. Otherwise, | |
| unknown number 1 will be associated with node 1. | |
| If IBC is 1 or 4, then there is an unknown value of U | |
| at node NSUB, which will be unknown NSUB or NSUB+1, | |
| depending on whether there was an unknown at node 0. | |
| Input, int NL. | |
| The number of basis functions used in a single | |
| subinterval. (NL-1) is the degree of the polynomials | |
| used. For this code, NL is fixed at 2, meaning that | |
| piecewise linear functions are used as the basis. | |
| Input, int NODE[NL*NSUB]. | |
| For each subinterval I: | |
| NODE[0+I*2] is the number of the left node, and | |
| NODE[1+I*2] is the number of the right node. | |
| Input, int NU. | |
| NU is the number of unknowns in the linear system. | |
| Depending on the value of IBC, there will be NSUB-1, | |
| NSUB, or NSUB+1 unknown values, which are the coefficients | |
| of basis functions. | |
| Input, int NQUAD. | |
| The number of quadrature points used in a subinterval. | |
| This code uses NQUAD = 1. | |
| Input, int NSUB. | |
| The number of subintervals into which the interval [XL,XR] is broken. | |
| Input, double UL. | |
| If IBC is 1 or 3, UL is the value that U is required | |
| to have at X = XL. | |
| If IBC is 2 or 4, UL is the value that U' is required | |
| to have at X = XL. | |
| Input, double UR. | |
| If IBC is 2 or 3, UR is the value that U is required | |
| to have at X = XR. | |
| If IBC is 1 or 4, UR is the value that U' is required | |
| to have at X = XR. | |
| Input, double XL. | |
| XL is the left endpoint of the interval over which the | |
| differential equation is being solved. | |
| Input, double XR. | |
| XR is the right endpoint of the interval over which the | |
| differential equation is being solved. | |
| */ | |
| { | |
| double aij; | |
| double he; | |
| int i; | |
| int ie; | |
| int ig; | |
| int il; | |
| int iq; | |
| int iu; | |
| int jg; | |
| int jl; | |
| int ju; | |
| double phii; | |
| double phiix; | |
| double phij; | |
| double phijx; | |
| double x; | |
| double xleft; | |
| double xquade; | |
| double xrite; | |
| /* | |
| Zero out the arrays that hold the coefficients of the matrix | |
| and the right hand side. | |
| */ | |
| for ( i = 0; i < nu; i++ ) | |
| { | |
| f[i] = 0.0; | |
| } | |
| for ( i = 0; i < nu; i++ ) | |
| { | |
| adiag[i] = 0.0; | |
| } | |
| for ( i = 0; i < nu; i++ ) | |
| { | |
| aleft[i] = 0.0; | |
| } | |
| for ( i = 0; i < nu; i++ ) | |
| { | |
| arite[i] = 0.0; | |
| } | |
| /* | |
| For interval number IE, | |
| */ | |
| for ( ie = 0; ie < nsub; ie++ ) | |
| { | |
| he = h[ie]; | |
| xleft = xn[node[0+ie*2]]; | |
| xrite = xn[node[1+ie*2]]; | |
| /* | |
| consider each quadrature point IQ, | |
| */ | |
| for ( iq = 0; iq < nquad; iq++ ) | |
| { | |
| xquade = xquad[ie]; | |
| /* | |
| and evaluate the integrals associated with the basis functions | |
| for the left, and for the right nodes. | |
| */ | |
| for ( il = 1; il <= nl; il++ ) | |
| { | |
| ig = node[il-1+ie*2]; | |
| iu = indx[ig] - 1; | |
| if ( 0 <= iu ) | |
| { | |
| phi ( il, xquade, &phii, &phiix, xleft, xrite ); | |
| f[iu] = f[iu] + he * ff ( xquade ) * phii; | |
| /* | |
| Take care of boundary nodes at which U' was specified. | |
| */ | |
| if ( ig == 0 ) | |
| { | |
| x = 0.0; | |
| f[iu] = f[iu] - pp ( x ) * ul; | |
| } | |
| else if ( ig == nsub ) | |
| { | |
| x = 1.0; | |
| f[iu] = f[iu] + pp ( x ) * ur; | |
| } | |
| /* | |
| Evaluate the integrals that take a product of the basis | |
| function times itself, or times the other basis function | |
| that is nonzero in this interval. | |
| */ | |
| for ( jl = 1; jl <= nl; jl++ ) | |
| { | |
| jg = node[jl-1+ie*2]; | |
| ju = indx[jg] - 1; | |
| phi ( jl, xquade, &phij, &phijx, xleft, xrite ); | |
| aij = he * ( pp ( xquade ) * phiix * phijx | |
| + qq ( xquade ) * phii * phij ); | |
| /* | |
| If there is no variable associated with the node, then it's | |
| a specified boundary value, so we multiply the coefficient | |
| times the specified boundary value and subtract it from the | |
| right hand side. | |
| */ | |
| if ( ju < 0 ) | |
| { | |
| if ( jg == 0 ) | |
| { | |
| f[iu] = f[iu] - aij * ul; | |
| } | |
| else if ( jg == nsub ) | |
| { | |
| f[iu] = f[iu] - aij * ur; | |
| } | |
| } | |
| /* | |
| Otherwise, we add the coefficient we've just computed to the | |
| diagonal, or left or right entries of row IU of the matrix. | |
| */ | |
| else | |
| { | |
| if ( iu == ju ) | |
| { | |
| adiag[iu] = adiag[iu] + aij; | |
| } | |
| else if ( ju < iu ) | |
| { | |
| aleft[iu] = aleft[iu] + aij; | |
| } | |
| else | |
| { | |
| arite[iu] = arite[iu] + aij; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| } | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ | |
| double ff ( double x ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| FF evaluates the right hand side function. | |
| Discussion: | |
| This routine evaluates the function F(X) in the differential equation. | |
| -d/dx (p du/dx) + q u = f | |
| at the point X. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, double X, the argument of the function. | |
| Output, double FF, the value of the function. | |
| */ | |
| { | |
| double value; | |
| value = 0.0; | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| void geometry ( double h[], int ibc, int indx[], int nl, int node[], int nsub, | |
| int *nu, double xl, double xn[], double xquad[], double xr ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| GEOMETRY sets up the geometry for the interval [XL,XR]. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Output, double H(NSUB) | |
| H(I) is the length of subinterval I. This code uses | |
| equal spacing for all the subintervals. | |
| Input, int IBC. | |
| IBC declares what the boundary conditions are. | |
| 1, at the left endpoint, U has the value UL, | |
| at the right endpoint, U' has the value UR. | |
| 2, at the left endpoint, U' has the value UL, | |
| at the right endpoint, U has the value UR. | |
| 3, at the left endpoint, U has the value UL, | |
| and at the right endpoint, U has the value UR. | |
| 4, at the left endpoint, U' has the value UL, | |
| at the right endpoint U' has the value UR. | |
| Output, int INDX[NSUB+1]. | |
| For a node I, INDX(I) is the index of the unknown | |
| associated with node I. | |
| If INDX(I) is equal to -1, then no unknown is associated | |
| with the node, because a boundary condition fixing the | |
| value of U has been applied at the node instead. | |
| Unknowns are numbered beginning with 1. | |
| If IBC is 2 or 4, then there is an unknown value of U | |
| at node 0, which will be unknown number 1. Otherwise, | |
| unknown number 1 will be associated with node 1. | |
| If IBC is 1 or 4, then there is an unknown value of U | |
| at node NSUB, which will be unknown NSUB or NSUB+1, | |
| depending on whether there was an unknown at node 0. | |
| Input, int NL. | |
| The number of basis functions used in a single | |
| subinterval. (NL-1) is the degree of the polynomials | |
| used. For this code, NL is fixed at 2, meaning that | |
| piecewise linear functions are used as the basis. | |
| Output, int NODE[NL*NSUB]. | |
| For each subinterval I: | |
| NODE[0+I*2] is the number of the left node, and | |
| NODE[1+I*2] is the number of the right node. | |
| Input, int NSUB. | |
| The number of subintervals into which the interval [XL,XR] is broken. | |
| Output, int *NU. | |
| NU is the number of unknowns in the linear system. | |
| Depending on the value of IBC, there will be NSUB-1, | |
| NSUB, or NSUB+1 unknown values, which are the coefficients | |
| of basis functions. | |
| Input, double XL. | |
| XL is the left endpoint of the interval over which the | |
| differential equation is being solved. | |
| Output, double XN(0:NSUB). | |
| XN(I) is the location of the I-th node. XN(0) is XL, | |
| and XN(NSUB) is XR. | |
| Output, double XQUAD(NSUB) | |
| XQUAD(I) is the location of the single quadrature point | |
| in interval I. | |
| Input, double XR. | |
| XR is the right endpoint of the interval over which the | |
| differential equation is being solved. | |
| */ | |
| { | |
| int i; | |
| /* | |
| Set the value of XN, the locations of the nodes. | |
| */ | |
| printf ( "\n" ); | |
| printf ( " Node Location\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i <= nsub; i++ ) | |
| { | |
| xn[i] = ( ( double ) ( nsub - i ) * xl | |
| + ( double ) i * xr ) | |
| / ( double ) ( nsub ); | |
| printf ( " %8d %14f \n", i, xn[i] ); | |
| } | |
| /* | |
| Set the lengths of each subinterval. | |
| */ | |
| printf ( "\n" ); | |
| printf ( "Subint Length\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i < nsub; i++ ) | |
| { | |
| h[i] = xn[i+1] - xn[i]; | |
| printf ( " %8d %14f\n", i+1, h[i] ); | |
| } | |
| /* | |
| Set the quadrature points, each of which is the midpoint | |
| of its subinterval. | |
| */ | |
| printf ( "\n" ); | |
| printf ( "Subint Quadrature point\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i < nsub; i++ ) | |
| { | |
| xquad[i] = 0.5 * ( xn[i] + xn[i+1] ); | |
| printf ( " %8d %14f\n", i+1, xquad[i] ); | |
| } | |
| /* | |
| Set the value of NODE, which records, for each interval, | |
| the node numbers at the left and right. | |
| */ | |
| printf ( "\n" ); | |
| printf ( "Subint Left Node Right Node\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i < nsub; i++ ) | |
| { | |
| node[0+i*2] = i; | |
| node[1+i*2] = i + 1; | |
| printf ( " %8d %8d %8d\n", i+1, node[0+i*2], node[1+i*2] ); | |
| } | |
| /* | |
| Starting with node 0, see if an unknown is associated with | |
| the node. If so, give it an index. | |
| */ | |
| *nu = 0; | |
| /* | |
| Handle first node. | |
| */ | |
| i = 0; | |
| if ( ibc == 1 || ibc == 3 ) | |
| { | |
| indx[i] = -1; | |
| } | |
| else | |
| { | |
| *nu = *nu + 1; | |
| indx[i] = *nu; | |
| } | |
| /* | |
| Handle nodes 1 through nsub-1 | |
| */ | |
| for ( i = 1; i < nsub; i++ ) | |
| { | |
| *nu = *nu + 1; | |
| indx[i] = *nu; | |
| } | |
| /* | |
| Handle the last node. | |
| /*/ | |
| i = nsub; | |
| if ( ibc == 2 || ibc == 3 ) | |
| { | |
| indx[i] = -1; | |
| } | |
| else | |
| { | |
| *nu = *nu + 1; | |
| indx[i] = *nu; | |
| } | |
| printf ( "\n" ); | |
| printf ( " Number of unknowns NU = %8d\n", *nu ); | |
| printf ( "\n" ); | |
| printf ( " Node Unknown\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i <= nsub; i++ ) | |
| { | |
| printf ( " %8d %8d\n", i, indx[i] ); | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ | |
| void init ( int *ibc, int *nquad, double *ul, double *ur, double *xl, | |
| double *xr ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| INIT assigns values to variables which define the problem. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Output, int *IBC. | |
| IBC declares what the boundary conditions are. | |
| 1, at the left endpoint, U has the value UL, | |
| at the right endpoint, U' has the value UR. | |
| 2, at the left endpoint, U' has the value UL, | |
| at the right endpoint, U has the value UR. | |
| 3, at the left endpoint, U has the value UL, | |
| and at the right endpoint, U has the value UR. | |
| 4, at the left endpoint, U' has the value UL, | |
| at the right endpoint U' has the value UR. | |
| Output, int *NQUAD. | |
| The number of quadrature points used in a subinterval. | |
| This code uses NQUAD = 1. | |
| Output, double *UL. | |
| If IBC is 1 or 3, UL is the value that U is required | |
| to have at X = XL. | |
| If IBC is 2 or 4, UL is the value that U' is required | |
| to have at X = XL. | |
| Output, double *UR. | |
| If IBC is 2 or 3, UR is the value that U is required | |
| to have at X = XR. | |
| If IBC is 1 or 4, UR is the value that U' is required | |
| to have at X = XR. | |
| Output, double *XL. | |
| XL is the left endpoint of the interval over which the | |
| differential equation is being solved. | |
| Output, double *XR. | |
| XR is the right endpoint of the interval over which the | |
| differential equation is being solved. | |
| */ | |
| { | |
| /* | |
| IBC declares what the boundary conditions are. | |
| */ | |
| *ibc = 1; | |
| /* | |
| NQUAD is the number of quadrature points per subinterval. | |
| The program as currently written cannot handle any value for | |
| NQUAD except 1. | |
| */ | |
| *nquad = 1; | |
| /* | |
| Set the values of U or U' at the endpoints. | |
| */ | |
| *ul = 0.0; | |
| *ur = 1.0; | |
| /* | |
| Define the location of the endpoints of the interval. | |
| */ | |
| *xl = 0.0; | |
| *xr = 1.0; | |
| /* | |
| Print out the values that have been set. | |
| */ | |
| printf ( "\n" ); | |
| printf ( " The equation is to be solved for\n" ); | |
| printf ( " X greater than XL = %f\n", *xl ); | |
| printf ( " and less than XR = %f\n", *xr ); | |
| printf ( "\n" ); | |
| printf ( " The boundary conditions are:\n" ); | |
| printf ( "\n" ); | |
| if ( *ibc == 1 || *ibc == 3 ) | |
| { | |
| printf ( " At X = XL, U = %f\n", *ul ); | |
| } | |
| else | |
| { | |
| printf ( " At X = XL, U' = %f\n", *ul ); | |
| } | |
| if ( *ibc == 2 || *ibc == 3 ) | |
| { | |
| printf ( " At X = XR, U = %f\n", *ur ); | |
| } | |
| else | |
| { | |
| printf ( " At X = XR, U' = %f\n", *ur ); | |
| } | |
| printf ( "\n" ); | |
| printf ( " Number of quadrature points per element is %d\n", *nquad ); | |
| return; | |
| } | |
| /******************************************************************************/ | |
| void output ( double f[], int ibc, int indx[], int nsub, int nu, double ul, | |
| double ur, double xn[] ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| OUTPUT prints out the computed solution. | |
| Discussion: | |
| We simply print out the solution vector F, except that, for | |
| certain boundary conditions, we are going to have to get the | |
| value of the solution at XL or XR by using the specified | |
| boundary value. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Input, double F(NU). | |
| ASSEMBLE stores into F the right hand side of the linear | |
| equations. | |
| SOLVE replaces those values of F by the solution of the | |
| linear equations. | |
| Input, int IBC. | |
| IBC declares what the boundary conditions are. | |
| 1, at the left endpoint, U has the value UL, | |
| at the right endpoint, U' has the value UR. | |
| 2, at the left endpoint, U' has the value UL, | |
| at the right endpoint, U has the value UR. | |
| 3, at the left endpoint, U has the value UL, | |
| and at the right endpoint, U has the value UR. | |
| 4, at the left endpoint, U' has the value UL, | |
| at the right endpoint U' has the value UR. | |
| Input, int INDX[NSUB+1]. | |
| For a node I, INDX(I) is the index of the unknown | |
| associated with node I. | |
| If INDX(I) is equal to -1, then no unknown is associated | |
| with the node, because a boundary condition fixing the | |
| value of U has been applied at the node instead. | |
| Unknowns are numbered beginning with 1. | |
| If IBC is 2 or 4, then there is an unknown value of U | |
| at node 0, which will be unknown number 1. Otherwise, | |
| unknown number 1 will be associated with node 1. | |
| If IBC is 1 or 4, then there is an unknown value of U | |
| at node NSUB, which will be unknown NSUB or NSUB+1, | |
| depending on whether there was an unknown at node 0. | |
| Input, int NSUB. | |
| The number of subintervals into which the interval [XL,XR] is broken. | |
| Input, int NU. | |
| NU is the number of unknowns in the linear system. | |
| Depending on the value of IBC, there will be NSUB-1, | |
| NSUB, or NSUB+1 unknown values, which are the coefficients | |
| of basis functions. | |
| Input, double UL. | |
| If IBC is 1 or 3, UL is the value that U is required | |
| to have at X = XL. | |
| If IBC is 2 or 4, UL is the value that U' is required | |
| to have at X = XL. | |
| Input, double UR. | |
| If IBC is 2 or 3, UR is the value that U is required | |
| to have at X = XR. | |
| If IBC is 1 or 4, UR is the value that U' is required | |
| to have at X = XR. | |
| Input, double XN(0:NSUB). | |
| XN(I) is the location of the I-th node. XN(0) is XL, | |
| and XN(NSUB) is XR. | |
| */ | |
| { | |
| int i; | |
| double u; | |
| printf ( "\n" ); | |
| printf ( " Computed solution coefficients:\n" ); | |
| printf ( "\n" ); | |
| printf ( " Node X(I) U(X(I))\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i <= nsub; i++ ) | |
| { | |
| /* | |
| If we're at the first node, check the boundary condition. | |
| */ | |
| if ( i == 0 ) | |
| { | |
| if ( ibc == 1 || ibc == 3 ) | |
| { | |
| u = ul; | |
| } | |
| else | |
| { | |
| u = f[indx[i]-1]; | |
| } | |
| } | |
| /* | |
| If we're at the last node, check the boundary condition. | |
| */ | |
| else if ( i == nsub ) | |
| { | |
| if ( ibc == 2 || ibc == 3 ) | |
| { | |
| u = ur; | |
| } | |
| else | |
| { | |
| u = f[indx[i]-1]; | |
| } | |
| } | |
| /* | |
| Any other node, we're sure the value is stored in F. | |
| */ | |
| else | |
| { | |
| u = f[indx[i]-1]; | |
| } | |
| printf ( " %8d %8f %14f\n", i, xn[i], u ); | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ | |
| void phi ( int il, double x, double *phii, double *phiix, double xleft, | |
| double xrite ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| PHI evaluates a linear basis function and its derivative. | |
| Discussion: | |
| The evaluation is done at a point X in an interval [XLEFT,XRITE]. | |
| In this interval, there are just two nonzero basis functions. | |
| The first basis function is a line which is 1 at the left | |
| endpoint and 0 at the right. The second basis function is 0 at | |
| the left endpoint and 1 at the right. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Input, int IL, the index of the basis function. | |
| 1, the function which is 1 at XLEFT and 0 at XRITE. | |
| 2, the function which is 0 at XLEFT and 1 at XRITE. | |
| Input, double X, the evaluation point. | |
| Output, double *PHII, *PHIIX, the value of the | |
| basis function and its derivative at X. | |
| Input, double XLEFT, XRITE, the left and right | |
| endpoints of the interval. | |
| */ | |
| { | |
| if ( xleft <= x && x <= xrite ) | |
| { | |
| if ( il == 1 ) | |
| { | |
| *phii = ( xrite - x ) / ( xrite - xleft ); | |
| *phiix = -1.0 / ( xrite - xleft ); | |
| } | |
| else | |
| { | |
| *phii = ( x - xleft ) / ( xrite - xleft ); | |
| *phiix = 1.0 / ( xrite - xleft ); | |
| } | |
| } | |
| /* | |
| If X is outside of the interval, just set everything to 0. | |
| */ | |
| else | |
| { | |
| *phii = 0.0; | |
| *phiix = 0.0; | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ | |
| double pp ( double x ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| PP evaluates the function P in the differential equation. | |
| Discussion: | |
| The function P appears in the differential equation as; | |
| - d/dx (p du/dx) + q u = f | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, double X, the argument of the function. | |
| Output, double PP, the value of the function. | |
| */ | |
| { | |
| double value; | |
| value = 1.0; | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| void prsys ( double adiag[], double aleft[], double arite[], double f[], | |
| int nu ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| PRSYS prints out the tridiagonal linear system. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameter: | |
| Input, double ADIAG(NU), the "diagonal" coefficients. That is, | |
| ADIAG(I) is the coefficient of the I-th unknown in the I-th equation. | |
| Input, double ALEFT(NU), the "left hand" coefficients. That is, ALEFT(I) | |
| is the coefficient of the (I-1)-th unknown in the I-th equation. | |
| There is no value in ALEFT(1), since the first equation | |
| does not refer to a "0-th" unknown. | |
| Input, double ARITE(NU). | |
| ARITE(I) is the "right hand" coefficient of the I-th | |
| equation in the linear system. ARITE(I) is the coefficient | |
| of the (I+1)-th unknown in the I-th equation. There is | |
| no value in ARITE(NU) because the NU-th equation does not | |
| refer to an "NU+1"-th unknown. | |
| Input, double F(NU). | |
| ASSEMBLE stores into F the right hand side of the linear | |
| equations. | |
| SOLVE replaces those values of F by the solution of the | |
| linear equations. | |
| Input, int NU. | |
| NU is the number of unknowns in the linear system. | |
| Depending on the value of IBC, there will be NSUB-1, | |
| NSUB, or NSUB+1 unknown values, which are the coefficients | |
| of basis functions. | |
| */ | |
| { | |
| int i; | |
| printf ( "\n" ); | |
| printf ( "Printout of tridiagonal linear system:\n" ); | |
| printf ( "\n" ); | |
| printf ( "Equation ALEFT ADIAG ARITE RHS\n" ); | |
| printf ( "\n" ); | |
| for ( i = 0; i < nu; i++ ) | |
| { | |
| printf ( " %8d %14f %14f %14f %14f\n", | |
| i + 1, aleft[i], adiag[i], arite[i], f[i] ); | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ | |
| double qq ( double x ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| QQ evaluates the function Q in the differential equation. | |
| Discussion: | |
| The function Q appears in the differential equation as: | |
| - d/dx (p du/dx) + q u = f | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| John Burkardt | |
| Parameters: | |
| Input, double X, the argument of the function. | |
| Output, double QQ, the value of the function. | |
| */ | |
| { | |
| double value; | |
| value = 0.0; | |
| return value; | |
| } | |
| /******************************************************************************/ | |
| void solve ( double adiag[], double aleft[], double arite[], double f[], | |
| int nu ) | |
| /******************************************************************************/ | |
| /* | |
| Purpose: | |
| SOLVE solves a tridiagonal matrix system of the form A*x = b. | |
| Licensing: | |
| This code is distributed under the GNU LGPL license. | |
| Modified: | |
| 29 May 2009 | |
| Author: | |
| C version by John Burkardt | |
| Parameters: | |
| Input/output, double ADIAG(NU), ALEFT(NU), ARITE(NU). | |
| On input, ADIAG, ALEFT, and ARITE contain the diagonal, | |
| left and right entries of the equations. | |
| On output, ADIAG and ARITE have been changed in order | |
| to compute the solution. | |
| Note that for the first equation, there is no ALEFT | |
| coefficient, and for the last, there is no ARITE. | |
| So there is no need to store a value in ALEFT(1), nor | |
| in ARITE(NU). | |
| Input/output, double F(NU). | |
| On input, F contains the right hand side of the linear | |
| system to be solved. | |
| On output, F contains the solution of the linear system. | |
| Input, int NU, the number of equations to be solved. | |
| */ | |
| { | |
| int i; | |
| /* | |
| Carry out Gauss elimination on the matrix, saving information | |
| needed for the backsolve. | |
| */ | |
| arite[0] = arite[0] / adiag[0]; | |
| for ( i = 1; i < nu - 1; i++ ) | |
| { | |
| adiag[i] = adiag[i] - aleft[i] * arite[i-1]; | |
| arite[i] = arite[i] / adiag[i]; | |
| } | |
| adiag[nu-1] = adiag[nu-1] - aleft[nu-1] * arite[nu-2]; | |
| /* | |
| Carry out the same elimination steps on F that were done to the | |
| matrix. | |
| */ | |
| f[0] = f[0] / adiag[0]; | |
| for ( i = 1; i < nu; i++ ) | |
| { | |
| f[i] = ( f[i] - aleft[i] * f[i-1] ) / adiag[i]; | |
| } | |
| /* | |
| And now carry out the steps of "back substitution". | |
| */ | |
| for ( i = nu - 2; 0 <= i; i-- ) | |
| { | |
| f[i] = f[i] - arite[i] * f[i+1]; | |
| } | |
| return; | |
| } | |
| /******************************************************************************/ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #! /usr/bin/env python3 | |
| # | |
| def barebones ( n, a, c, f, x ): | |
| #*****************************************************************************80 | |
| # | |
| ## barebones ??? | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/barebones.py | |
| # | |
| import numpy as np | |
| import scipy.linalg as la | |
| quad_num = 2 | |
| abscissa = np.array ( [ -0.577350269189625764509148780502, \ | |
| +0.577350269189625764509148780502 ] ) | |
| weight = np.array ( [ 1.0, 1.0 ] ) | |
| A = np.zeros ( [ n, n ] ) | |
| b = np.zeros ( n ) | |
| e_num = n - 1 | |
| for e in range ( 0, e_num ): | |
| l = e | |
| xl = x[l] | |
| r = e + 1 | |
| xr = x[r] | |
| for q in range ( 0, quad_num ): | |
| xq = ( ( 1.0 - abscissa[q] ) * xl \ | |
| + ( 1.0 + abscissa[q] ) * xr ) \ | |
| / 2.0 | |
| wq = weight[q] * ( xr - xl ) / 2.0 | |
| vl = ( xr - xq ) / ( xr - xl ) | |
| vlp = - 1.0 / ( xr - xl ) | |
| vr = ( xq - xl ) / ( xr - xl ) | |
| vrp = +1.0 / ( xr - xl ) | |
| axq = a ( xq ) | |
| cxq = c ( xq ) | |
| fxq = f ( xq ) | |
| A[l,l] = A[l,l] + wq * ( vlp * axq * vlp + vl * cxq * vl ) | |
| A[l,r] = A[l,r] + wq * ( vlp * axq * vrp + vl * cxq * vr ) | |
| b[l] = b[l] + wq * ( vl * fxq ) | |
| A[r,l] = A[r,l] + wq * ( vrp * axq * vlp + vr * cxq * vl ) | |
| A[r,r] = A[r,r] + wq * ( vrp * axq * vrp + vr * cxq * vr ) | |
| b[r] = b[r] + wq * ( vr * fxq ) | |
| for j in range ( 0, n ): | |
| A[0,j] = 0.0 | |
| A[0,0] = 1.0 | |
| b[0] = 0.0 | |
| for j in range ( 0, n ): | |
| A[n-1,j] = 0.0 | |
| A[n-1,n-1] = 1.0 | |
| b[n-1] = 0.0 | |
| u = la.solve ( A, b ) | |
| return u | |
| def barebones_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## barebones_test tests barebones. | |
| # | |
| import matplotlib.pyplot as plt | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'BAREBONES_TEST' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A(X) = 1.0' ) | |
| print ( ' C(X) = 1.0' ) | |
| print ( ' F(X) = X' ) | |
| print ( ' U(X) = X - SINH(X) / SINH(1)' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| x_lo = 0.0 | |
| x_hi = 1.0 | |
| x = np.linspace ( x_lo, x_hi, n ) | |
| u = barebones ( n, a, c, f, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact ) | |
| e2 = l2_error_linear ( n, x, u, exact ) | |
| h1s = h1s_error_linear ( n, x, u, exactp ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| fig = plt.figure ( ) | |
| plt.plot ( x, u, 'bo-' ) | |
| plt.xlabel ( '<---X--->' ) | |
| plt.ylabel ( '<---U(X)--->' ) | |
| plt.title ( 'BAREBONES Solution' ) | |
| plt.savefig ( 'barebones.png' ) | |
| plt.show ( block = False ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'BAREBONES_TEST' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a ( x ): | |
| #*****************************************************************************80 | |
| # | |
| value = 1.0 | |
| return value | |
| def c ( x ): | |
| #*****************************************************************************80 | |
| # | |
| value = 1.0 | |
| return value | |
| def exact ( x ): | |
| #*****************************************************************************80 | |
| # | |
| import numpy as np | |
| value = x - np.sinh ( x ) / np.sinh ( 1.0 ) | |
| return value | |
| def exactp ( x ): | |
| #*****************************************************************************80 | |
| # | |
| import numpy as np | |
| value = 1.0 - np.cosh ( x ) / np.sinh ( 1.0 ) | |
| return value | |
| def f ( x ): | |
| #*****************************************************************************80 | |
| # | |
| value = x | |
| return value | |
| def fem1d_bvp_linear ( n, a, c, f, x ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR solves a two point boundary value problem. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear.py | |
| # | |
| # Discussion: | |
| # | |
| # The program uses the finite element method, with piecewise linear basis | |
| # functions to solve a boundary value problem in one dimension. | |
| # | |
| # The problem is defined on the region 0 <= x <= 1. | |
| # | |
| # The following differential equation is imposed between 0 and 1: | |
| # | |
| # - d/dx a(x) du/dx + c(x) * u(x) = f(x) | |
| # | |
| # where a(x), c(x), and f(x) are given functions. | |
| # | |
| # At the boundaries, the following conditions are applied: | |
| # | |
| # u(0.0) = 0.0 | |
| # u(1.0) = 0.0 | |
| # | |
| # A set of N equally spaced nodes is defined on this | |
| # interval, with 0 = X(1) < X(2) < ... < X(N) = 1.0. | |
| # | |
| # At each node I, we associate a piecewise linear basis function V(I,X), | |
| # which is 0 at all nodes except node I. This implies that V(I,X) is | |
| # everywhere 0 except that | |
| # | |
| # for X(I-1) <= X <= X(I): | |
| # | |
| # V(I,X) = ( X - X(I-1) ) / ( X(I) - X(I-1) ) | |
| # | |
| # for X(I) <= X <= X(I+1): | |
| # | |
| # V(I,X) = ( X(I+1) - X ) / ( X(I+1) - X(I) ) | |
| # | |
| # We now assume that the solution U(X) can be written as a linear | |
| # sum of these basis functions: | |
| # | |
| # U(X) = sum ( 1 <= J <= N ) U(J) * V(J,X) | |
| # | |
| # where U(X) on the left is the function of X, but on the right, | |
| # is meant to indicate the coefficients of the basis functions. | |
| # | |
| # To determine the coefficient U(J), we multiply the original | |
| # differential equation by the basis function V(J,X), and use | |
| # integration by parts, to arrive at the I-th finite element equation: | |
| # | |
| # Integral A(X) * U'(X) * V'(I,X) + C(X) * U(X) * V(I,X) dx | |
| # = Integral F(X) * V(I,X) dx | |
| # | |
| # We note that the functions U(X) and U'(X) can be replaced by | |
| # the finite element form involving the linear sum of basis functions, | |
| # but we also note that the resulting integrand will only be nonzero | |
| # for terms where J = I - 1, I, or I + 1. | |
| # | |
| # By writing this equation for basis functions I = 2 through N - 1, | |
| # and using the boundary conditions, we have N linear equations | |
| # for the N unknown coefficients U(1) through U(N), which can | |
| # be easily solved. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # Input, integer N, the number of nodes. | |
| # | |
| # Input, function A ( X ), evaluates a(x); | |
| # | |
| # Input, function C ( X ), evaluates c(x); | |
| # | |
| # Input, function F ( X ), evaluates f(x); | |
| # | |
| # Input, real X(N), the mesh points. | |
| # | |
| # Output, real U(N), the finite element coefficients, which are also | |
| # the value of the computed solution at the mesh points. | |
| # | |
| import numpy as np | |
| import scipy.linalg as la | |
| # | |
| # Define a 2 point Gauss-Legendre quadrature rule on [-1,+1]. | |
| # | |
| quad_num = 2 | |
| abscissa = np.array ( [ -0.577350269189625764509148780502, \ | |
| +0.577350269189625764509148780502 ] ) | |
| weight = np.array ( [ 1.0, 1.0 ] ) | |
| # | |
| # Make room for the matrix A and right hand side b. | |
| # | |
| A = np.zeros ( [ n, n ] ) | |
| b = np.zeros ( n ) | |
| # | |
| # We assemble the finite element matrix by looking at each element, | |
| # that is, each interval [ X(L), X(R) ]. | |
| # | |
| # There are only two nonzero piecewise linear basis functions in this | |
| # element, and we can call them VL and VR. So the only integrals we | |
| # need to compute involve products of: | |
| # | |
| # VL * VL VR * VL F * VL | |
| # VR * VL VR * VR F * VR | |
| # | |
| # and | |
| # | |
| # VL' * VL' VR' * VL' | |
| # VR' * VL' VR' * VR' | |
| # | |
| e_num = n - 1 | |
| for e in range ( 0, e_num ): | |
| l = e | |
| xl = x[l] | |
| r = e + 1 | |
| xr = x[r] | |
| for q in range ( 0, quad_num ): | |
| xq = ( ( 1.0 - abscissa[q] ) * xl \ | |
| + ( 1.0 + abscissa[q] ) * xr ) \ | |
| / 2.0 | |
| wq = weight[q] * ( xr - xl ) / 2.0 | |
| vl = ( xr - xq ) / ( xr - xl ) | |
| vlp = - 1.0 / ( xr - xl ) | |
| vr = ( xq - xl ) / ( xr - xl ) | |
| vrp = +1.0 / ( xr - xl ) | |
| axq = a ( xq ) | |
| cxq = c ( xq ) | |
| fxq = f ( xq ) | |
| A[l,l] = A[l,l] + wq * ( vlp * axq * vlp + vl * cxq * vl ) | |
| A[l,r] = A[l,r] + wq * ( vlp * axq * vrp + vl * cxq * vr ) | |
| b[l] = b[l] + wq * ( vl * fxq ) | |
| A[r,l] = A[r,l] + wq * ( vrp * axq * vlp + vr * cxq * vl ) | |
| A[r,r] = A[r,r] + wq * ( vrp * axq * vrp + vr * cxq * vr ) | |
| b[r] = b[r] + wq * ( vr * fxq ) | |
| # | |
| # Equation 0 is the left boundary condition, U(0) = 0.0; | |
| # | |
| for j in range ( 0, n ): | |
| A[0,j] = 0.0 | |
| A[0,0] = 1.0 | |
| b[0] = 0.0 | |
| # | |
| # We can keep the matrix symmetric by using the boundary condition | |
| # to eliminate U(0) from all equations but #0. | |
| # | |
| for i in range ( 1, n ): | |
| b[i] = b[i] - A[i,0] * b[0] | |
| A[i,0] = 0.0 | |
| # | |
| # Equation N-1 is the right boundary condition, U(N-1) = 0.0; | |
| # | |
| for j in range ( 0, n ): | |
| A[n-1,j] = 0.0 | |
| A[n-1,n-1] = 1.0 | |
| b[n-1] = 0.0 | |
| # | |
| # We can keep the matrix symmetric by using the boundary condition | |
| # to eliminate U(N-1) from all equations but #N-1. | |
| # | |
| for i in range ( 0, n - 1 ): | |
| b[i] = b[i] - A[i,n-1] * b[n-1] | |
| A[i,n-1] = 0.0 | |
| # | |
| # Solve the linear system for the finite element coefficients U. | |
| # | |
| u = la.solve ( A, b ) | |
| return u | |
| def fem1d_bvp_linear_test00 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST00 tests FEM1D_BVP_LINEAR. | |
| # | |
| # 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) | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/m_src/fem1d_bvp_linear/fem1d_bvp_linear_test00.py | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import matplotlib.pyplot as plt | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST00' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A(X) = 1.0' ) | |
| print ( ' C(X) = 1.0' ) | |
| print ( ' F(X) = X' ) | |
| print ( ' U(X) = X - SINH(X) / SINH(1)' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_lo = 0.0 | |
| x_hi = 1.0 | |
| x = np.linspace ( x_lo, x_hi, n ) | |
| u = fem1d_bvp_linear ( n, a00, c00, f00, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact00 ( x[i] ) | |
| # | |
| # Print a table. | |
| # | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| # | |
| # Compute error norms. | |
| # | |
| e1 = l1_error ( n, x, u, exact00 ) | |
| e2 = l2_error_linear ( n, x, u, exact00 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp00 ) | |
| mx = max_error_linear ( n, x, u, exact00 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Plot the computed solution. | |
| # | |
| fig = plt.figure ( ) | |
| plt.plot ( x, u, 'bo-' ) | |
| plt.xlabel ( '<---X--->' ) | |
| plt.ylabel ( '<---U(X)--->' ) | |
| plt.title ( 'FEM1D_BVP_LINEAR_TEST00 Solution' ) | |
| plt.savefig ( 'fem1d_bvp_linear_test00.png' ) | |
| plt.show ( block = False ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST00' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a00 ( x ): | |
| value = 1.0 | |
| return value | |
| def c00 ( x ): | |
| value = 1.0 | |
| return value | |
| def exact00 ( x ): | |
| from math import sinh | |
| value = x - sinh ( x ) / sinh ( 1.0 ) | |
| return value | |
| def exactp00 ( x ): | |
| from math import cosh | |
| from math import sinh | |
| value = 1.0 - cosh ( x ) / sinh ( 1.0 ) | |
| return value | |
| def f00 ( x ): | |
| value = x | |
| return value | |
| def fem1d_bvp_linear_test01 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST01 carries out test case #1. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test01.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A1, C1, F1, EXACT1, EXACT_UX1. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST01' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A1(X) = 1.0' ) | |
| print ( ' C1(X) = 0.0' ) | |
| print ( ' F1(X) = X * ( X + 3 ) * exp ( X )' ) | |
| print ( ' U1(X) = X * ( 1 - X ) * exp ( X )' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a1, c1, f1, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact1 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact1 ) | |
| e2 = l2_error_linear ( n, x, u, exact1 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp1 ) | |
| mx = max_error_linear ( n, x, u, exact1 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST01' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a1 ( x ): | |
| value = 1.0 | |
| return value | |
| def c1 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact1 ( x ): | |
| from math import exp | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| return value | |
| def exactp1 ( x ): | |
| from math import exp | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| return value | |
| def f1 ( x ): | |
| from math import exp | |
| value = x * ( x + 3.0 ) * exp ( x ) | |
| return value | |
| def fem1d_bvp_linear_test02 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST02 carries out test case #2. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test02.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A2, C2, F2, EXACT2, EXACT_UX2. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST02' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A2(X) = 1.0' ) | |
| print ( ' C2(X) = 2.0' ) | |
| print ( ' F2(X) = X * ( 5 - X ) * exp ( X )' ) | |
| print ( ' U2(X) = X * ( 1 - X ) * exp ( X )' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a2, c2, f2, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact2 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact2 ) | |
| e2 = l2_error_linear ( n, x, u, exact2 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp2 ) | |
| mx = max_error_linear ( n, x, u, exact2 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST02' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a2 ( x ): | |
| value = 1.0 | |
| return value | |
| def c2 ( x ): | |
| value = 2.0 | |
| return value | |
| def exact2 ( x ): | |
| from math import exp | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| return value | |
| def exactp2 ( x ): | |
| from math import exp | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| return value | |
| def f2 ( x ): | |
| from math import exp | |
| value = x * ( 5.0 - x ) * exp ( x ) | |
| return value | |
| def fem1d_bvp_linear_test03 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST03 carries out test case #3. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test03.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A3, C3, F3, EXACT3, EXACT_UX3. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST03' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A3(X) = 1.0' ) | |
| print ( ' C3(X) = 2.0 * X' ) | |
| print ( ' F3(X) = - X * ( 2 * X * X - 3 * X - 3 ) * exp ( X )' ) | |
| print ( ' U3(X) = X * ( 1 - X ) * exp ( X )' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a3, c3, f3, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact3 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact3 ) | |
| e2 = l2_error_linear ( n, x, u, exact3 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp3 ) | |
| mx = max_error_linear ( n, x, u, exact3 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST03' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a3 ( x ): | |
| value = 1.0 | |
| return value | |
| def c3 ( x ): | |
| value = 2.0 * x | |
| return value | |
| def exact3( x ): | |
| from math import exp | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| return value | |
| def exactp3 ( x ): | |
| from math import exp | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| return value | |
| def f3 ( x ): | |
| from math import exp | |
| value = - x * ( 2.0 * x * x - 3.0 * x - 3.0 ) * exp ( x ) | |
| return value | |
| def fem1d_bvp_linear_test04 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST04 carries out test case #4. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test04.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A4, C4, F4, EXACT4, EXACT_UX4. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST04' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A4(X) = 1.0 + X * X' ) | |
| print ( ' C4(X) = 0.0' ) | |
| print ( ' F4(X) = ( X + 3 X^2 + 5 X^3 + X^4 ) * exp ( X )' ) | |
| print ( ' U4(X) = X * ( 1 - X ) * exp ( X )' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a4, c4, f4, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact4 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact4 ) | |
| e2 = l2_error_linear ( n, x, u, exact4 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp4 ) | |
| mx = max_error_linear ( n, x, u, exact4 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST04' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a4 ( x ): | |
| value = 1.0 + x * x | |
| return value | |
| def c4 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact4 ( x ): | |
| from math import exp | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| return value | |
| def exactp4 ( x ): | |
| from math import exp | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| return value | |
| def f4 ( x ): | |
| from math import exp | |
| value = ( x + 3.0 * x * x + 5.0 * x ** 3 + x ** 4 ) * exp ( x ) | |
| return value | |
| def fem1d_bvp_linear_test05 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST05 carries out test case #5. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test05.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A5, C5, F5, EXACT5, EXACT_UX5. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST05' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A5(X) = 1.0 + X * X for X <= 1/3' ) | |
| print ( ' = 7/9 + X for 1/3 < X' ) | |
| print ( ' C5(X) = 0.0' ) | |
| print ( ' F5(X) = ( X + 3 X^2 + 5 X^3 + X^4 ) * exp ( X )' ) | |
| print ( ' for X <= 1/3' ) | |
| print ( ' = ( - 1 + 10/3 X + 43/9 X^2 + X^3 ) .* exp ( X )' ) | |
| print ( ' for 1/3 <= X' ) | |
| print ( ' U5(X) = X * ( 1 - X ) * exp ( X )' ) | |
| print ( '' ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a5, c5, f5, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact5 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact5 ) | |
| e2 = l2_error_linear ( n, x, u, exact5 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp5 ) | |
| mx = max_error_linear ( n, x, u, exact5 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST05' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a5 ( x ): | |
| if ( x <= 1.0 / 3.0 ): | |
| value = 1.0 + x * x | |
| else: | |
| value = x + 7.0 / 9.0 | |
| return value | |
| def c5 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact5 ( x ): | |
| from math import exp | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| return value | |
| def exactp5 ( x ): | |
| from math import exp | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| return value | |
| def f5 ( x ): | |
| from math import exp | |
| if ( x <= 1.0 / 3.0 ): | |
| value = ( x + 3.0 * x ** 2 + 5.0 * x ** 3 + x ** 4 ) * exp ( x ) | |
| else: | |
| value = ( - 1.0 + ( 10.0 / 3.0 ) * x \ | |
| + ( 43.0 / 9.0 ) * x ** 2 + x ** 3 ) * exp ( x ) | |
| return value | |
| def fem1d_bvp_linear_test06 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST06 does an error analysis. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test06.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A6, C6, F6, EXACT6, EXACT_UX6. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST06' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A6(X) = 1.0' ) | |
| print ( ' C6(X) = 0.0' ) | |
| print ( ' F6(X) = pi*pi*sin(pi*X)' ) | |
| print ( ' U6(X) = sin(pi*x)' ) | |
| print ( '' ) | |
| print ( ' Compute L2 norm and seminorm of error for various N.' ) | |
| print ( '' ) | |
| print ( ' N L1 error L2 error Seminorm error Maxnorm error' ) | |
| print ( '' ) | |
| n = 11 | |
| for i in range ( 0, 5 ): | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| 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, exactp6 ) | |
| mx = max_error_linear ( n, x, u, exact6 ) | |
| print ( ' %4d %14g %14g %14g %14g' % ( n, e1, e2, h1s, mx ) ) | |
| n = 2 * ( n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST06' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a6 ( x ): | |
| value = 1.0 | |
| return value | |
| def c6 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact6 ( x ): | |
| import numpy as np | |
| value = np.sin ( np.pi * x ) | |
| return value | |
| def exactp6 ( x ): | |
| import numpy as np | |
| value = np.pi * np.cos ( np.pi * x ) | |
| return value | |
| def f6 ( x ): | |
| import numpy as np | |
| value = np.pi ** 2 * np.sin ( np.pi * x ) | |
| return value | |
| def fem1d_bvp_linear_test07 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST07 does an error analysis. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test07.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A7, C7, F7, EXACT7, EXACT_UX7. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # 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. | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST07' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Becker/Carey/Oden example.' ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( '' ) | |
| print ( ' Compute L2 norm and seminorm of error for various N.' ) | |
| print ( '' ) | |
| print ( ' N L1 error L2 error Seminorm error Maxnorm error' ) | |
| print ( '' ) | |
| n = 11 | |
| for i in range ( 0, 5 ): | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| 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, exactp7 ) | |
| mx = max_error_linear ( n, x, u, exact7 ) | |
| print ( ' %4d %14g %14g %14g %14g' % ( n, e1, e2, h1s, mx ) ) | |
| n = 2 * ( n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST07' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a7 ( x ): | |
| alpha = 30.0 | |
| x0 = 1.0 / 3.0 | |
| value = 1.0 / alpha + alpha * ( x - x0 ) ** 2 | |
| return value | |
| def c7 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact7 ( x ): | |
| from math import atan | |
| alpha = 30.0 | |
| x0 = 1.0 / 3.0 | |
| value = ( 1.0 - x ) * ( atan ( alpha * ( x - x0 ) ) + atan ( alpha * x0 ) ) | |
| return value | |
| def exactp7 ( x ): | |
| from math import atan | |
| alpha = 30.0 | |
| x0 = 1.0 / 3.0 | |
| value = - atan ( alpha * ( x - x0 ) ) - atan ( alpha * x0 ) \ | |
| + ( 1.0 - x ) * alpha / ( 1.0 + alpha * alpha * ( x - x0 ) ** 2 ) | |
| return value | |
| def f7 ( x ): | |
| from math import atan | |
| 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 | |
| def fem1d_bvp_linear_test08 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST08 carries out test case #8. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test08.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A8, C8, F8, EXACT8, EXACT_UX8. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST08' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A8(X) = 1.0' ) | |
| print ( ' C8(X) = 0.0' ) | |
| print ( ' F8(X) = X * ( X + 3 ) * exp ( X ), X <= 2/3' ) | |
| print ( ' = 2 * exp ( 2/3), 2/3 < X' ) | |
| print ( ' U8(X) = X * ( 1 - X ) * exp ( X ), X <= 2/3' ) | |
| print ( ' = X * ( 1 - X ) * exp ( 2/3 ), 2/3 < X' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a8, c8, f8, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact8 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact8 ) | |
| e2 = l2_error_linear ( n, x, u, exact8 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp8 ) | |
| mx = max_error_linear ( n, x, u, exact8 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST08' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a8 ( x ): | |
| value = 1.0 | |
| return value | |
| def c8 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact8 ( x ): | |
| from math import exp | |
| 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 | |
| def exactp8 ( x ): | |
| from math import exp | |
| 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 | |
| def f8 ( x ): | |
| from math import exp | |
| if ( x <= 2.0 / 3.0 ): | |
| value = x * ( x + 3.0 ) * exp ( x ) | |
| else: | |
| value = 2.0 * exp ( 2.0 / 3.0 ) | |
| return value | |
| def fem1d_bvp_linear_test09 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST09 carries out test case #9. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test09.py | |
| # | |
| # Discussion: | |
| # | |
| # Use A9, C9, F9, EXACT9, EXACT_UX9. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Reference: | |
| # | |
| # Dianne O'Leary, | |
| # Scientific Computing with Case Studies, | |
| # SIAM, 2008, | |
| # ISBN13: 978-0-898716-66-5, | |
| # LC: QA401.O44. | |
| # | |
| import numpy as np | |
| import platform | |
| n = 11 | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST09' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A9(X) = 1.0' ) | |
| print ( ' C9(X) = 0.0' ) | |
| print ( ' F9(X) = X * ( X + 3 ) * exp ( X ), X <= 2/3' ) | |
| print ( ' = 2 * exp ( 2/3), 2/3 < X' ) | |
| print ( ' U9(X) = X * ( 1 - X ) * exp ( X ), X <= 2/3' ) | |
| print ( ' = X * ( 1 - X ), 2/3 < X' ) | |
| print ( '' ) | |
| print ( ' Number of nodes = %d' % ( n ) ) | |
| # | |
| # Geometry definitions. | |
| # | |
| x_first = 0.0 | |
| x_last = 1.0 | |
| x = np.linspace ( x_first, x_last, n ) | |
| u = fem1d_bvp_linear ( n, a9, c9, f9, x ) | |
| g = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| g[i] = exact9 ( x[i] ) | |
| print ( '' ) | |
| print ( ' I X U Uexact Error' ) | |
| print ( '' ) | |
| for i in range ( 0, n ): | |
| print ( ' %4d %8f %8f %8f %8e' \ | |
| % ( i, x[i], u[i], g[i], abs ( u[i] - g[i] ) ) ) | |
| e1 = l1_error ( n, x, u, exact9 ) | |
| e2 = l2_error_linear ( n, x, u, exact9 ) | |
| h1s = h1s_error_linear ( n, x, u, exactp9 ) | |
| mx = max_error_linear ( n, x, u, exact9 ) | |
| print ( '' ) | |
| print ( ' l1 norm of error = %g' % ( e1 ) ) | |
| print ( ' L2 norm of error = %g' % ( e2 ) ) | |
| print ( ' Seminorm of error = %g' % ( h1s ) ) | |
| print ( ' Max norm of error = %g' % ( mx ) ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST09' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a9 ( x ): | |
| value = 1.0 | |
| return value | |
| def c9 ( x ): | |
| value = 0.0 | |
| return value | |
| def exact9 ( x ): | |
| from math import exp | |
| if ( x <= 2.0 / 3.0 ): | |
| value = x * ( 1.0 - x ) * exp ( x ) | |
| else: | |
| value = x * ( 1.0 - x ) | |
| return value | |
| def exactp9 ( x ): | |
| from math import exp | |
| if ( x <= 2.0 / 3.0 ): | |
| value = ( 1.0 - x - x * x ) * exp ( x ) | |
| else: | |
| value = ( 1.0 - 2.0 * x ) | |
| return value | |
| def f9 ( x ): | |
| from math import exp | |
| if ( x <= 2.0 / 3.0 ): | |
| value = x * ( x + 3.0 ) * exp ( x ) | |
| else: | |
| value = 2.0 | |
| return value | |
| def fem1d_bvp_linear_test10 ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_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 GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 16 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. | |
| # | |
| import matplotlib.pyplot as plt | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST10' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Solve -( A(x) U\'(x) )\' + C(x) U(x) = F(x)' ) | |
| print ( ' for 0 < x < 1, with U(0) = U(1) = 0.' ) | |
| print ( ' A(X) = 1.0' ) | |
| print ( ' C(X) = 1.0' ) | |
| print ( ' F(X) = X' ) | |
| print ( ' U(X) = X - SINH(X) / SINH(1)' ) | |
| print ( '' ) | |
| print ( ' log(E) E L2error H1error Maxerror' ) | |
| print ( '' ) | |
| e_log_max = 6 | |
| ne_plot = np.zeros ( e_log_max + 1 ) | |
| h_plot = np.zeros ( e_log_max + 1 ) | |
| l2_plot = np.zeros ( e_log_max + 1 ) | |
| h1_plot = np.zeros ( e_log_max + 1 ) | |
| mx_plot = np.zeros ( e_log_max + 1 ) | |
| for e_log in range ( 0, e_log_max + 1 ): | |
| ne = 2 ** e_log | |
| n = ne + 1 | |
| x_lo = 0.0 | |
| x_hi = 1.0 | |
| x = np.linspace ( x_lo, x_hi, n ) | |
| u = fem1d_bvp_linear ( n, a10, c10, f10, x ) | |
| ne_plot[e_log] = ne | |
| h_plot[e_log] = ( x_hi - x_lo ) / float ( ne ) | |
| l2_plot[e_log] = l2_error_linear ( n, x, u, exact10 ) | |
| h1_plot[e_log] = h1s_error_linear ( n, x, u, exactp10 ) | |
| mx_plot[e_log] = max_error_linear ( n, x, u, exact10 ) | |
| print ( ' %4d %4d %14g %14g %14g' \ | |
| % ( e_log, ne, l2_plot[e_log], h1_plot[e_log], mx_plot[e_log] ) ) | |
| print ( '' ) | |
| print ( ' log(E1) E1 / E2 L2rate H1rate Maxrate' ) | |
| print ( '' ) | |
| for e_log in range ( 0, e_log_max ): | |
| ne1 = ne_plot[e_log] | |
| ne2 = ne_plot[e_log+1] | |
| r = float ( ne2 ) / float ( ne1 ) | |
| l2 = l2_plot[e_log] / l2_plot[e_log+1] | |
| l2 = np.log ( l2 ) / np.log ( r ) | |
| h1 = h1_plot[e_log] / h1_plot[e_log+1] | |
| h1 = np.log ( h1 ) / np.log ( r ) | |
| mx = mx_plot[e_log] / mx_plot[e_log+1] | |
| mx = np.log ( mx ) / np.log ( r ) | |
| print ( ' %4d %4d/%4d %14g %14g %14g' \ | |
| % ( e_log, ne1, ne2, l2, h1, mx ) ) | |
| # | |
| # Plot the L2 error as a function of NE. | |
| # | |
| fig = plt.figure ( ) | |
| plt.plot ( ne_plot, l2_plot, 'bo-' ) | |
| plt.xlabel ( '<---NE--->' ) | |
| plt.ylabel ( '<---L2(error)--->' ) | |
| plt.title ( 'L2 error as function of number of elements' ) | |
| plt.xscale ( 'log' ) | |
| plt.yscale ( 'log' ) | |
| plt.grid ( True ) | |
| plt.axis ( 'equal' ) | |
| plt.savefig ( 'l2error_test10.png' ) | |
| plt.show ( block = False ) | |
| # | |
| # Plot the max error as a function of NE. | |
| # | |
| fig = plt.figure ( ) | |
| plt.plot ( ne_plot, mx_plot, 'bo-' ) | |
| plt.xlabel ( '<---NE--->' ) | |
| plt.ylabel ( '<---Max(error)--->' ) | |
| plt.title ( 'Max error as function of number of elements' ) | |
| plt.xscale ( 'log' ) | |
| plt.yscale ( 'log' ) | |
| plt.grid ( True ) | |
| plt.axis ( 'equal' ) | |
| plt.savefig ( 'maxerror_test10.png' ) | |
| plt.show ( block = False ) | |
| # | |
| # Plot the H1 error as a function of NE. | |
| # | |
| fig = plt.figure ( ) | |
| plt.plot ( ne_plot, h1_plot, 'bo-' ) | |
| plt.xlabel ( '<---NE--->' ) | |
| plt.ylabel ( '<---H1(error)--->' ) | |
| plt.title ( 'H1 error as function of number of elements' ) | |
| plt.xscale ( 'log' ) | |
| plt.yscale ( 'log' ) | |
| plt.grid ( True ) | |
| plt.axis ( 'equal' ) | |
| plt.savefig ( 'h1error_test10.png' ) | |
| plt.show ( block = False ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST10' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def a10 ( x ): | |
| value = 1.0 | |
| return value | |
| def c10 ( x ): | |
| value = 1.0 | |
| return value | |
| def exact10 ( x ): | |
| from math import sinh | |
| value = x - sinh ( x ) / sinh ( 1.0 ) | |
| return value | |
| def exactp10 ( x ): | |
| from math import cosh | |
| from math import sinh | |
| value = 1.0 - cosh ( x ) / sinh ( 1.0 ) | |
| return value | |
| def f10 ( x ): | |
| value = x | |
| return value | |
| def fem1d_bvp_linear_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## FEM1D_BVP_LINEAR_TEST tests the FEM1D_BVP_LINEAR library. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/fem1d_bvp_linear_test.py | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 07 July 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| import platform | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' Test the FEM1D_BVP_LINEAR library.' ) | |
| # | |
| # Utility functions. | |
| # | |
| h1s_error_linear_test ( ) | |
| l1_error_test ( ) | |
| l2_error_linear_test ( ) | |
| max_error_linear_test ( ) | |
| # | |
| # Library functions. | |
| # | |
| fem1d_bvp_linear_test00 ( ) | |
| fem1d_bvp_linear_test01 ( ) | |
| fem1d_bvp_linear_test02 ( ) | |
| fem1d_bvp_linear_test03 ( ) | |
| fem1d_bvp_linear_test04 ( ) | |
| fem1d_bvp_linear_test05 ( ) | |
| fem1d_bvp_linear_test06 ( ) | |
| fem1d_bvp_linear_test07 ( ) | |
| fem1d_bvp_linear_test08 ( ) | |
| fem1d_bvp_linear_test09 ( ) | |
| fem1d_bvp_linear_test10 ( ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'FEM1D_BVP_LINEAR_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| def h1s_error_linear ( n, x, u, exact_ux ): | |
| #*****************************************************************************80 | |
| # | |
| ## H1S_ERROR_LINEAR: seminorm error of a finite element solution. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/h1s_error_linear.py | |
| # | |
| # Discussion: | |
| # | |
| # We assume the finite element method has been used, over an interval [A,B] | |
| # involving N nodes, with piecewise linear elements used for the basis. | |
| # The finite element solution U(x) has been computed, and a formula for the | |
| # exact derivative V'(x) is known. | |
| # | |
| # This function estimates the H1 seminorm of the error: | |
| # | |
| # H1S = sqrt ( integral ( A <= x <= B ) ( U'(x) - V'(x) )^2 dx | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 12 July 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # Input, integer N, the number of nodes. | |
| # | |
| # Input, real X(N), the mesh points. | |
| # | |
| # Input, real U(N), the finite element coefficients. | |
| # | |
| # Input, function EQ = EXACT_UX ( X ), returns the value of the exact | |
| # derivative at the point X. | |
| # | |
| # Output, real H1S, the estimated seminorm of the error. | |
| # | |
| import numpy as np | |
| h1s = 0.0 | |
| # | |
| # Define a 2 point Gauss-Legendre quadrature rule on [-1,+1]. | |
| # | |
| quad_num = 2 | |
| abscissa = np.array ( [ -0.577350269189625764509148780502, \ | |
| +0.577350269189625764509148780502 ] ) | |
| weight = np.array ( [ 1.0, 1.0 ] ) | |
| # | |
| # Integrate over each interval. | |
| # | |
| e_num = n - 1 | |
| for e in range ( 0, e_num ): | |
| l = e | |
| xl = x[l] | |
| ul = u[l] | |
| r = e + 1 | |
| xr = x[r] | |
| ur = u[r] | |
| for q in range ( 0, quad_num ): | |
| xq = ( ( 1.0 - abscissa[q] ) * xl \ | |
| + ( 1.0 + abscissa[q] ) * xr ) \ | |
| / 2.0 | |
| wq = weight[q] * ( xr - xl ) / 2.0 | |
| # | |
| # The piecewise linear derivative is a constant in the interval. | |
| # | |
| uxq = ( ur - ul ) / ( xr - xl ) | |
| exq = exact_ux ( xq ) | |
| h1s = h1s + wq * ( uxq - exq ) ** 2 | |
| h1s = np.sqrt ( h1s ) | |
| return h1s | |
| def h1s_error_linear_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## H1S_ERROR_LINEAR_TEST tests H1S_ERROR_LINEAR. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'H1S_ERROR_LINEAR_TEST:' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' H1S_ERROR_LINEAR computes the H1 seminorm approximation error' ) | |
| print ( ' between the exact derivative of a function and the derivative' ) | |
| print ( ' of a piecewise linear approximation to the function,' ) | |
| print ( ' associated with n mesh points x().' ) | |
| print ( '' ) | |
| print ( ' N H1S_Error' ) | |
| print ( '' ) | |
| x_n = 3 | |
| for test in range ( 0, 8 ): | |
| x_lo = 0.0 | |
| x_hi = np.pi | |
| x = np.linspace ( x_lo, x_hi, x_n ) | |
| # | |
| # U is an approximation to sin(x). | |
| # | |
| u = np.zeros ( x_n ) | |
| for i in range ( 0, x_n ): | |
| u[i] = np.sin ( x[i] ) | |
| # | |
| # Compare the derivative of the piecewise interpolant of U | |
| # to the actual derivative, cos(x). | |
| # | |
| e1 = h1s_error_linear ( x_n, x, u, np.cos ) | |
| print ( ' %2d %g' % ( x_n, e1 ) ) | |
| x_n = 2 * ( x_n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'H1S_ERROR_LINEAR_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def l1_error ( n, x, u, exact ): | |
| #*****************************************************************************80 | |
| # | |
| ## L1_ERROR estimates the l1 error norm of a finite element solution. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/l1_error.py | |
| # | |
| # Discussion: | |
| # | |
| # We assume the finite element method has been used, over an interval [A,B] | |
| # involving N nodes. | |
| # | |
| # The coefficients U(1:N) have been computed, and a formula for the | |
| # exact solution is known. | |
| # | |
| # This function estimates the little l1 norm of the error: | |
| # L1_NORM = sum ( 1 <= I <= N ) abs ( U(i) - EXACT(X(i)) ) | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # Input, integer N, the number of nodes. | |
| # | |
| # Input, real X(N), the mesh points. | |
| # | |
| # Input, real U(N), the finite element coefficients. | |
| # | |
| # Input, function EQ = EXACT ( X ), returns the value of the exact | |
| # solution at the point X. | |
| # | |
| # Output, real E1, the little l1 norm of the error. | |
| # | |
| e1 = 0.0 | |
| for i in range ( 0, n ): | |
| e1 = e1 + abs ( u[i] - exact ( x[i] ) ) | |
| e1 = e1 / n | |
| return e1 | |
| def l1_error_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## L1_ERROR_TEST tests L1_ERROR. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'L1_ERROR_TEST:' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' L1_ERROR computes the little l1 approximation error between' ) | |
| print ( ' a function exact(x) and a vector of n values u() at points x().' ) | |
| print ( '' ) | |
| print ( ' N L1_Error' ) | |
| print ( '' ) | |
| x_n = 3 | |
| for test in range ( 0, 6 ): | |
| x_lo = 0.0 | |
| x_hi = np.pi | |
| x = np.linspace ( x_lo, x_hi, x_n ) | |
| # | |
| # U is an approximation to sin(x). | |
| # | |
| u = np.zeros ( x_n ) | |
| for i in range ( 0, x_n ): | |
| u[i] = x[i] - x[i] ** 3 / 6.0 | |
| e1 = l1_error ( x_n, x, u, np.sin ) | |
| print ( ' %2d %g' % ( x_n, e1 ) ) | |
| x_n = 2 * ( x_n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'L1_ERROR_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def l2_error_linear ( n, x, u, exact ): | |
| #*****************************************************************************80 | |
| # | |
| ## L2_ERROR_LINEAR estimates the L2 error norm of a finite element solution. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/l2_error_linear.py | |
| # | |
| # Discussion: | |
| # | |
| # We assume the finite element method has been used, over an interval [A,B] | |
| # involving N nodes, with piecewise linear elements used for the basis. | |
| # The coefficients U(1:N) have been computed, and a formula for the | |
| # exact solution is known. | |
| # | |
| # This function estimates the L2 norm of the error: | |
| # | |
| # L2_NORM = Integral ( A <= X <= B ) ( U(X) - EXACT(X) )^2 dX | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 12 July 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # Input, integer N, the number of nodes. | |
| # | |
| # Input, real X(N), the mesh points. | |
| # | |
| # Input, real U(N), the finite element coefficients. | |
| # | |
| # Input, function EQ = EXACT ( X ), returns the value of the exact | |
| # solution at the point X. | |
| # | |
| # Output, real E2, the estimated L2 norm of the error. | |
| # | |
| import numpy as np | |
| e2 = 0.0 | |
| # | |
| # Define a 2 point Gauss-Legendre quadrature rule on [-1,+1]. | |
| # | |
| quad_num = 2 | |
| abscissa = np.array ( [ -0.577350269189625764509148780502, \ | |
| +0.577350269189625764509148780502 ] ) | |
| weight = np.array ( [ 1.0, 1.0 ] ) | |
| # | |
| # Integrate over each interval. | |
| # | |
| e_num = n - 1 | |
| for e in range ( 0, e_num ): | |
| l = e | |
| xl = x[l] | |
| ul = u[l] | |
| r = e + 1 | |
| xr = x[r] | |
| ur = u[r] | |
| for q in range ( 0, quad_num ): | |
| xq = ( ( 1.0 - abscissa[q] ) * xl \ | |
| + ( 1.0 + abscissa[q] ) * xr ) \ | |
| / 2.0 | |
| wq = weight[q] * ( xr - xl ) / 2.0 | |
| # | |
| # Use the fact that U is a linear combination of piecewise linears. | |
| # | |
| uq = ( ( xr - xq ) * ul \ | |
| + ( xq - xl ) * ur ) \ | |
| / ( xr - xl ) | |
| eq = exact ( xq ) | |
| e2 = e2 + wq * ( uq - eq ) ** 2 | |
| e2 = np.sqrt ( e2 ) | |
| return e2 | |
| def l2_error_linear_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## L2_ERROR_LINEAR_TEST tests L2_ERROR_LINEAR. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 09 January 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'L2_ERROR_LINEAR_TEST:' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' L2_ERROR_LINEAR computes the L2 approximation error between' ) | |
| print ( ' a function exact(x) and a piecewise linear function u()' ) | |
| print ( ' associated with n mesh points x().' ) | |
| print ( '' ) | |
| print ( ' N L2_Error' ) | |
| print ( '' ) | |
| x_n = 3 | |
| for test in range ( 0, 6 ): | |
| x_lo = 0.0 | |
| x_hi = np.pi | |
| x = np.linspace ( x_lo, x_hi, x_n ) | |
| # | |
| # U is an approximation to sin(x). | |
| # | |
| u = np.zeros ( x_n ) | |
| for i in range ( 0, x_n ): | |
| u[i] = np.sin ( x[i] ) | |
| e1 = l2_error_linear ( x_n, x, u, np.sin ) | |
| print ( ' %2d %g' % ( x_n, e1 ) ) | |
| x_n = 2 * ( x_n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'L2_ERROR_LINEAR_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def max_error_linear ( n, x, u, exact ): | |
| #*****************************************************************************80 | |
| # | |
| ## MAX_ERROR_LINEAR estimates the max error norm of a finite element solution. | |
| # | |
| # Location: | |
| # | |
| # http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/max_error_linear.py | |
| # | |
| # Discussion: | |
| # | |
| # We assume the finite element method has been used, over an interval [A,B] | |
| # involving N nodes, with piecewise linear elements used for the basis. | |
| # The coefficients U(1:N) have been computed, and a formula for the | |
| # exact solution is known. | |
| # | |
| # This function estimates the max norm of the error: | |
| # | |
| # MAX_NORM = Integral ( A <= X <= B ) max ( abs ( U(X) - EXACT(X) ) ) dX | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 07 July 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # Input, integer N, the number of nodes. | |
| # | |
| # Input, real X(N), the mesh points. | |
| # | |
| # Input, real U(N), the finite element coefficients. | |
| # | |
| # Input, function EQ = EXACT ( X ), returns the value of the exact | |
| # solution at the point X. | |
| # | |
| # Output, real VALUE, the estimated max norm of the error. | |
| # | |
| import numpy as np | |
| quad_num = 8 | |
| value = 0.0 | |
| # | |
| # Examine QUAD_NUM points in each element, including left node but not right. | |
| # | |
| e_num = n - 1 | |
| for e in range ( 0, e_num ): | |
| l = e | |
| xl = x[l] | |
| ul = u[l] | |
| r = e + 1 | |
| xr = x[r] | |
| ur = u[r] | |
| for q in range ( 0, quad_num ): | |
| xq = ( float ( quad_num - q ) * xl \ | |
| + float ( q ) * xr ) \ | |
| / float ( quad_num ) | |
| # | |
| # Use the fact that U is a linear combination of piecewise linears. | |
| # | |
| uq = ( ( xr - xq ) * ul \ | |
| + ( xq - xl ) * ur ) \ | |
| / ( xr - xl ) | |
| eq = exact ( xq ) | |
| value = max ( value, abs ( uq - eq ) ) | |
| # | |
| # For completeness, check last node. | |
| # | |
| xq = x[n-1] | |
| uq = u[n-1] | |
| eq = exact ( xq ) | |
| value = max ( value, abs ( uq - eq ) ) | |
| # | |
| # Integral approximation requires multiplication by interval length. | |
| # | |
| value = value * ( x[n-1] - x[0] ) | |
| return value | |
| def max_error_linear_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## MAX_ERROR_LINEAR_TEST tests MAX_ERROR_LINEAR. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 07 July 2015 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import numpy as np | |
| import platform | |
| print ( '' ) | |
| print ( 'MAX_ERROR_LINEAR_TEST:' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' MAX_ERROR_LINEAR computes the maximum absolute approximation error' ) | |
| print ( ' between a function exact(x) and a piecewise linear function u()' ) | |
| print ( ' associated with n mesh points x().' ) | |
| print ( '' ) | |
| print ( ' N Max_Error' ) | |
| print ( '' ) | |
| n = 3 | |
| for test in range ( 0, 5 ): | |
| x_lo = 0.0 | |
| x_hi = np.pi | |
| x = np.linspace ( x_lo, x_hi, n ) | |
| # | |
| # U is an approximation to sin(x). | |
| # | |
| u = np.zeros ( n ) | |
| for i in range ( 0, n ): | |
| u[i] = np.sin ( x[i] ) | |
| e1 = max_error_linear ( n, x, u, np.sin ) | |
| print ( ' %2d %g' % ( n, e1 ) ) | |
| n = 2 * ( n - 1 ) + 1 | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'MAX_ERROR_LINEAR_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| def timestamp ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## TIMESTAMP prints the date as a timestamp. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 06 April 2013 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import time | |
| t = time.time ( ) | |
| print ( time.ctime ( t ) ) | |
| return None | |
| def timestamp_test ( ): | |
| #*****************************************************************************80 | |
| # | |
| ## TIMESTAMP_TEST tests TIMESTAMP. | |
| # | |
| # Licensing: | |
| # | |
| # This code is distributed under the GNU LGPL license. | |
| # | |
| # Modified: | |
| # | |
| # 03 December 2014 | |
| # | |
| # Author: | |
| # | |
| # John Burkardt | |
| # | |
| # Parameters: | |
| # | |
| # None | |
| # | |
| import platform | |
| print ( '' ) | |
| print ( 'TIMESTAMP_TEST:' ) | |
| print ( ' Python version: %s' % ( platform.python_version ( ) ) ) | |
| print ( ' TIMESTAMP prints a timestamp of the current date and time.' ) | |
| print ( '' ) | |
| timestamp ( ) | |
| # | |
| # Terminate. | |
| # | |
| print ( '' ) | |
| print ( 'TIMESTAMP_TEST:' ) | |
| print ( ' Normal end of execution.' ) | |
| return | |
| if ( __name__ == '__main__' ): | |
| timestamp ( ) | |
| fem1d_bvp_linear_test ( ) | |
| timestamp ( ) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include <stdio.h> | |
| #include <stdlib.h> | |
| #include <math.h> | |
| #include <time.h> | |
| #define LEARNING_RATE 0.1 | |
| #define MAX_ITERATION 100 | |
| float randomFloat() | |
| { | |
| return (float)rand() / (float)RAND_MAX; | |
| } | |
| int calculateOutput(float weights[], float x, float y) | |
| { | |
| float sum = x * weights[0] + y * weights[1] + weights[2]; | |
| return (sum >= 0) ? 1 : -1; | |
| } | |
| int main(int argc, char *argv[]) | |
| { | |
| srand(time(NULL)); | |
| float x[208], y[208], weights[3], localError, globalError; | |
| int outputs[208], patternCount, i, p, iteration, output; | |
| FILE *fp; | |
| if ((fp = fopen("test1.txt", "r")) == NULL) { | |
| printf("Cannot open file.\n"); | |
| exit(1); | |
| } | |
| i = 0; | |
| while (fscanf(fp, "%f %f %d", &x[i], &y[i], &outputs[i]) != EOF) { | |
| if (outputs[i] == 0) { | |
| outputs[i] = -1; | |
| } | |
| i++; | |
| } | |
| patternCount = i; | |
| weights[0] = randomFloat(); | |
| weights[1] = randomFloat(); | |
| weights[2] = randomFloat(); | |
| iteration = 0; | |
| do { | |
| iteration++; | |
| globalError = 0; | |
| for (p = 0; p < patternCount; p++) { | |
| output = calculateOutput(weights, x[p], y[p]); | |
| localError = outputs[p] - output; | |
| weights[0] += LEARNING_RATE * localError * x[p]; | |
| weights[1] += LEARNING_RATE * localError * y[p]; | |
| weights[2] += LEARNING_RATE * localError; | |
| globalError += (localError*localError); | |
| } | |
| /* Root Mean Squared Error */ | |
| printf("Iteration %d : RMSE = %.4f\n", iteration, | |
| sqrt(globalError/patternCount)); | |
| } while (globalError != 0 && iteration<=MAX_ITERATION); | |
| printf("\nDecision boundary (line) equation: %.2f*x + %.2f*y + %.2f = 0\n", | |
| weights[0], weights[1], weights[2]); | |
| return 0; | |
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #https://github.com/hayderkharrufa/snake_pathfinding_ai | |
| # Dimensions | |
| WIDTH = 612 # Width of game surface | |
| HEIGHT = 612 # Height of game surface | |
| ROWS = 17 | |
| SQUARE_SIZE = WIDTH // ROWS | |
| GAP_SIZE = 2 # Gap between adjacent squares | |
| # Colors | |
| SURFACE_CLR = (15, 15, 15) | |
| GRID_CLR = (20, 20, 20) | |
| SNAKE_CLR = (50, 255, 50) | |
| APPLE_CLR = (255, 255, 0) | |
| HEAD_CLR = (0, 150, 0) | |
| VIRTUAL_SNAKE_CLR = (255, 0, 0) | |
| # Game Settings | |
| FPS = 30 # Frames per second | |
| INITIAL_SNAKE_LENGTH = 3 | |
| WAIT_SECONDS_AFTER_WIN = 15 # If snake wins the game, wait for this amount of seconds before restarting | |
| MAX_MOVES_WITHOUT_EATING = ROWS * ROWS * ROWS * 2 # Snake will die after this amount of moves without eating apple | |
| SNAKE_MAX_LENGTH = ROWS * ROWS - INITIAL_SNAKE_LENGTH # Max number of apples snake can eat | |
| # Variables used in BFS algorithm | |
| GRID = [[i, j] for i in range(ROWS) for j in range(ROWS)] | |
| # Helper functions | |
| def get_neighbors(position): | |
| neighbors = [[position[0] + 1, position[1]], | |
| [position[0] - 1, position[1]], | |
| [position[0], position[1] + 1], | |
| [position[0], position[1] - 1]] | |
| in_grid_neighbors = [] | |
| for pos in neighbors: | |
| if pos in GRID: | |
| in_grid_neighbors.append(pos) | |
| return in_grid_neighbors | |
| def distance(pos1, pos2): | |
| x1, x2 = pos1[0], pos2[0] | |
| y1, y2 = pos1[1], pos2[1] | |
| return abs(x2 - x1) + abs(y2 - y1) | |
| # Each position is a tuple because python doesn't allow hashing lists | |
| ADJACENCY_DICT = {tuple(pos): get_neighbors(pos) for pos in GRID} | |
| # settings 可以直接放入 brython 執行 | |
| #print(distance([0,0], [10,10])) | |
| from copy import deepcopy | |
| from random import randrange | |
| class Square: | |
| def __init__(self, pos, surface, is_apple=False): | |
| self.pos = pos | |
| self.surface = surface | |
| self.is_apple = is_apple | |
| self.is_tail = False | |
| self.dir = [-1, 0] # [x, y] Direction | |
| if self.is_apple: | |
| self.dir = [0, 0] | |
| def draw(self, clr=SNAKE_CLR): | |
| x, y = self.pos[0], self.pos[1] | |
| ss, gs = SQUARE_SIZE, GAP_SIZE | |
| if self.dir == [-1, 0]: | |
| if self.is_tail: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss - 2*gs))''') | |
| else: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss, ss - 2*gs))''') | |
| if self.dir == [1, 0]: | |
| if self.is_tail: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss - 2*gs))''') | |
| else: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss - gs, y * ss + gs, ss, ss - 2*gs))''') | |
| if self.dir == [0, 1]: | |
| if self.is_tail: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss - 2*gs))''') | |
| else: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss - gs, ss - 2*gs, ss))''') | |
| if self.dir == [0, -1]: | |
| if self.is_tail: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss - 2*gs))''') | |
| else: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss))''') | |
| if self.is_apple: | |
| print('''pygame.draw.rect(self.surface, clr, (x * ss + gs, y * ss + gs, ss - 2*gs, ss - 2*gs))''') | |
| def move(self, direction): | |
| self.dir = direction | |
| self.pos[0] += self.dir[0] | |
| self.pos[1] += self.dir[1] | |
| def hitting_wall(self): | |
| if (self.pos[0] <= -1) or (self.pos[0] >= ROWS) or (self.pos[1] <= -1) or (self.pos[1] >= ROWS): | |
| return True | |
| else: | |
| return False | |
| class Snake: | |
| def __init__(self, surface): | |
| self.surface = surface | |
| self.is_dead = False | |
| self.squares_start_pos = [[ROWS // 2 + i, ROWS // 2] for i in range(INITIAL_SNAKE_LENGTH)] | |
| self.turns = {} | |
| self.dir = [-1, 0] | |
| self.score = 0 | |
| self.moves_without_eating = 0 | |
| self.apple = Square([randrange(ROWS), randrange(ROWS)], self.surface, is_apple=True) | |
| self.squares = [] | |
| for pos in self.squares_start_pos: | |
| self.squares.append(Square(pos, self.surface)) | |
| self.head = self.squares[0] | |
| self.tail = self.squares[-1] | |
| self.tail.is_tail = True | |
| self.path = [] | |
| self.is_virtual_snake = False | |
| self.total_moves = 0 | |
| self.won_game = False | |
| def draw(self): | |
| self.apple.draw(APPLE_CLR) | |
| self.head.draw(HEAD_CLR) | |
| for sqr in self.squares[1:]: | |
| if self.is_virtual_snake: | |
| sqr.draw(VIRTUAL_SNAKE_CLR) | |
| else: | |
| sqr.draw() | |
| def set_direction(self, direction): | |
| if direction == 'left': | |
| if not self.dir == [1, 0]: | |
| self.dir = [-1, 0] | |
| self.turns[self.head.pos[0], self.head.pos[1]] = self.dir | |
| if direction == "right": | |
| if not self.dir == [-1, 0]: | |
| self.dir = [1, 0] | |
| self.turns[self.head.pos[0], self.head.pos[1]] = self.dir | |
| if direction == "up": | |
| if not self.dir == [0, 1]: | |
| self.dir = [0, -1] | |
| self.turns[self.head.pos[0], self.head.pos[1]] = self.dir | |
| if direction == "down": | |
| if not self.dir == [0, -1]: | |
| self.dir = [0, 1] | |
| self.turns[self.head.pos[0], self.head.pos[1]] = self.dir | |
| def handle_events(self): | |
| for event in pygame.event.get(): | |
| if event.type == pygame.QUIT: | |
| pygame.quit() | |
| # Set snake direction using keyboard | |
| keys = pygame.key.get_pressed() | |
| for _ in keys: | |
| if keys[pygame.K_LEFT]: | |
| self.set_direction('left') | |
| elif keys[pygame.K_RIGHT]: | |
| self.set_direction('right') | |
| elif keys[pygame.K_UP]: | |
| self.set_direction('up') | |
| elif keys[pygame.K_DOWN]: | |
| self.set_direction('down') | |
| def move(self): | |
| for j, sqr in enumerate(self.squares): | |
| p = (sqr.pos[0], sqr.pos[1]) | |
| if p in self.turns: | |
| turn = self.turns[p] | |
| sqr.move([turn[0], turn[1]]) | |
| if j == len(self.squares) - 1: | |
| self.turns.pop(p) | |
| else: | |
| sqr.move(sqr.dir) | |
| self.moves_without_eating += 1 | |
| def add_square(self): | |
| self.squares[-1].is_tail = False | |
| tail = self.squares[-1] # Tail before adding new square | |
| direction = tail.dir | |
| if direction == [1, 0]: | |
| self.squares.append(Square([tail.pos[0] - 1, tail.pos[1]], self.surface)) | |
| if direction == [-1, 0]: | |
| self.squares.append(Square([tail.pos[0] + 1, tail.pos[1]], self.surface)) | |
| if direction == [0, 1]: | |
| self.squares.append(Square([tail.pos[0], tail.pos[1] - 1], self.surface)) | |
| if direction == [0, -1]: | |
| self.squares.append(Square([tail.pos[0], tail.pos[1] + 1], self.surface)) | |
| self.squares[-1].dir = direction | |
| self.squares[-1].is_tail = True # Tail after adding new square | |
| def reset(self): | |
| self.__init__(self.surface) | |
| def hitting_self(self): | |
| for sqr in self.squares[1:]: | |
| if sqr.pos == self.head.pos: | |
| return True | |
| def generate_apple(self): | |
| self.apple = Square([randrange(ROWS), randrange(ROWS)], self.surface, is_apple=True) | |
| if not self.is_position_free(self.apple.pos): | |
| self.generate_apple() | |
| def eating_apple(self): | |
| if self.head.pos == self.apple.pos and not self.is_virtual_snake and not self.won_game: | |
| self.generate_apple() | |
| self.moves_without_eating = 0 | |
| self.score += 1 | |
| return True | |
| def go_to(self, position): # Set head direction to target position | |
| if self.head.pos[0] - 1 == position[0]: | |
| self.set_direction('left') | |
| if self.head.pos[0] + 1 == position[0]: | |
| self.set_direction('right') | |
| if self.head.pos[1] - 1 == position[1]: | |
| self.set_direction('up') | |
| if self.head.pos[1] + 1 == position[1]: | |
| self.set_direction('down') | |
| def is_position_free(self, position): | |
| if position[0] >= ROWS or position[0] < 0 or position[1] >= ROWS or position[1] < 0: | |
| return False | |
| for sqr in self.squares: | |
| if sqr.pos == position: | |
| return False | |
| return True | |
| # Breadth First Search Algorithm | |
| def bfs(self, s, e): # Find shortest path between (start_position, end_position) | |
| q = [s] # Queue | |
| visited = {tuple(pos): False for pos in GRID} | |
| visited[s] = True | |
| # Prev is used to find the parent node of each node to create a feasible path | |
| prev = {tuple(pos): None for pos in GRID} | |
| while q: # While queue is not empty | |
| node = q.pop(0) | |
| neighbors = ADJACENCY_DICT[node] | |
| for next_node in neighbors: | |
| if self.is_position_free(next_node) and not visited[tuple(next_node)]: | |
| q.append(tuple(next_node)) | |
| visited[tuple(next_node)] = True | |
| prev[tuple(next_node)] = node | |
| path = list() | |
| p_node = e # Starting from end node, we will find the parent node of each node | |
| start_node_found = False | |
| while not start_node_found: | |
| if prev[p_node] is None: | |
| return [] | |
| p_node = prev[p_node] | |
| if p_node == s: | |
| path.append(e) | |
| return path | |
| path.insert(0, p_node) | |
| return [] # Path not available | |
| def create_virtual_snake(self): # Creates a copy of snake (same size, same position, etc..) | |
| v_snake = Snake(self.surface) | |
| for i in range(len(self.squares) - len(v_snake.squares)): | |
| v_snake.add_square() | |
| for i, sqr in enumerate(v_snake.squares): | |
| sqr.pos = deepcopy(self.squares[i].pos) | |
| sqr.dir = deepcopy(self.squares[i].dir) | |
| v_snake.dir = deepcopy(self.dir) | |
| v_snake.turns = deepcopy(self.turns) | |
| v_snake.apple.pos = deepcopy(self.apple.pos) | |
| v_snake.apple.is_apple = True | |
| v_snake.is_virtual_snake = True | |
| return v_snake | |
| def get_path_to_tail(self): | |
| tail_pos = deepcopy(self.squares[-1].pos) | |
| self.squares.pop(-1) | |
| path = self.bfs(tuple(self.head.pos), tuple(tail_pos)) | |
| self.add_square() | |
| return path | |
| def get_available_neighbors(self, pos): | |
| valid_neighbors = [] | |
| neighbors = get_neighbors(tuple(pos)) | |
| for n in neighbors: | |
| if self.is_position_free(n) and self.apple.pos != n: | |
| valid_neighbors.append(tuple(n)) | |
| return valid_neighbors | |
| def longest_path_to_tail(self): | |
| neighbors = self.get_available_neighbors(self.head.pos) | |
| path = [] | |
| if neighbors: | |
| dis = -9999 | |
| for n in neighbors: | |
| if distance(n, self.squares[-1].pos) > dis: | |
| v_snake = self.create_virtual_snake() | |
| v_snake.go_to(n) | |
| v_snake.move() | |
| if v_snake.eating_apple(): | |
| v_snake.add_square() | |
| if v_snake.get_path_to_tail(): | |
| path.append(n) | |
| dis = distance(n, self.squares[-1].pos) | |
| if path: | |
| return [path[-1]] | |
| def any_safe_move(self): | |
| neighbors = self.get_available_neighbors(self.head.pos) | |
| path = [] | |
| if neighbors: | |
| path.append(neighbors[randrange(len(neighbors))]) | |
| v_snake = self.create_virtual_snake() | |
| for move in path: | |
| v_snake.go_to(move) | |
| v_snake.move() | |
| if v_snake.get_path_to_tail(): | |
| return path | |
| else: | |
| return self.get_path_to_tail() | |
| def set_path(self): | |
| # If there is only 1 apple left for snake to win and it's adjacent to head | |
| if self.score == SNAKE_MAX_LENGTH - 1 and self.apple.pos in get_neighbors(self.head.pos): | |
| winning_path = [tuple(self.apple.pos)] | |
| print('Snake is about to win..') | |
| return winning_path | |
| v_snake = self.create_virtual_snake() | |
| # Let the virtual snake check if path to apple is available | |
| path_1 = v_snake.bfs(tuple(v_snake.head.pos), tuple(v_snake.apple.pos)) | |
| # This will be the path to virtual snake tail after it follows path_1 | |
| path_2 = [] | |
| if path_1: | |
| for pos in path_1: | |
| v_snake.go_to(pos) | |
| v_snake.move() | |
| v_snake.add_square() # Because it will eat an apple | |
| path_2 = v_snake.get_path_to_tail() | |
| # v_snake.draw() | |
| if path_2: # If there is a path between v_snake and it's tail | |
| return path_1 # Choose BFS path to apple (Fastest and shortest path) | |
| # If path_1 or path_2 not available, test these 3 conditions: | |
| # 1- Make sure that the longest path to tail is available | |
| # 2- If score is even, choose longest_path_to_tail() to follow the tail, if odd use any_safe_move() | |
| # 3- Change the follow tail method if the snake gets stuck in a loop | |
| if self.longest_path_to_tail() and\ | |
| self.score % 2 == 0 and\ | |
| self.moves_without_eating < MAX_MOVES_WITHOUT_EATING / 2: | |
| # Choose longest path to tail | |
| return self.longest_path_to_tail() | |
| # Play any possible safe move and make sure path to tail is available | |
| if self.any_safe_move(): | |
| return self.any_safe_move() | |
| # If path to tail is available | |
| if self.get_path_to_tail(): | |
| # Choose shortest path to tail | |
| return self.get_path_to_tail() | |
| # Snake couldn't find a path and will probably die | |
| print('No available path, snake in danger!') | |
| def update(self): | |
| self.handle_events() | |
| self.path = self.set_path() | |
| if self.path: | |
| self.go_to(self.path[0]) | |
| self.draw() | |
| self.move() | |
| if self.score == ROWS * ROWS - INITIAL_SNAKE_LENGTH: # If snake wins the game | |
| self.won_game = True | |
| print("Snake won the game after {} moves" | |
| .format(self.total_moves)) | |
| print('''pygame.time.wait(1000 * WAIT_SECONDS_AFTER_WIN)''') | |
| return 1 | |
| self.total_moves += 1 | |
| if self.hitting_self() or self.head.hitting_wall(): | |
| print("Snake is dead, trying again..") | |
| self.is_dead = True | |
| self.reset() | |
| if self.moves_without_eating == MAX_MOVES_WITHOUT_EATING: | |
| self.is_dead = True | |
| print("Snake got stuck, trying again..") | |
| self.reset() | |
| if self.eating_apple(): | |
| self.add_square() | |
| print("debug") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 1.456044 -7.058824 1 | |
| -4.478022 -2.072829 1 | |
| -7.664835 -6.890756 1 | |
| -5.137363 2.352941 0 | |
| -8.818681 3.025210 1 | |
| 3.653846 -2.969188 0 | |
| -5.686813 -2.128852 1 | |
| -5.961538 -2.577031 1 | |
| 5.961538 3.865546 0 | |
| -6.620879 0.168067 1 | |
| -8.159341 -5.434174 1 | |
| 5.027473 4.873950 0 | |
| 2.994505 -7.563025 1 | |
| -3.818681 -8.347339 1 | |
| 5.192308 0.336134 0 | |
| -8.653846 -2.184874 1 | |
| 4.148352 2.352941 0 | |
| 0.521978 -3.193277 0 | |
| -4.807692 -8.459384 1 | |
| -3.818681 -3.361345 1 | |
| -7.115385 1.120448 1 | |
| -8.763736 -4.761905 1 | |
| 5.192308 -8.067227 1 | |
| 3.598901 -2.296919 0 | |
| -5.576923 2.801120 0 | |
| -8.489011 -6.666667 1 | |
| 1.456044 6.498599 0 | |
| 0.796703 -5.154062 1 | |
| -2.994505 -8.907563 1 | |
| -2.280220 -5.378151 1 | |
| -3.708791 -8.907563 1 | |
| -5.247253 -4.649860 1 | |
| -0.851648 6.946779 0 | |
| -3.269231 -5.154062 1 | |
| -7.005495 -1.904762 1 | |
| -5.192308 -6.330532 1 | |
| -1.840659 -8.907563 1 | |
| -4.807692 -2.633053 1 | |
| -3.543956 3.865546 0 | |
| 4.917582 -3.809524 0 | |
| 1.401099 4.817927 0 | |
| 1.785714 -3.361345 0 | |
| -1.895604 2.689076 0 | |
| 7.390110 2.857143 0 | |
| -1.565934 -6.274510 1 | |
| 1.016484 -8.067227 1 | |
| -4.587912 -4.257703 1 | |
| -2.939560 -6.610644 1 | |
| -4.203297 6.946779 0 | |
| 3.269231 -8.851541 1 | |
| -3.049451 -5.714286 1 | |
| 4.587912 1.120448 0 | |
| -2.445055 5.658263 0 | |
| 3.818681 -9.019608 1 | |
| -2.060440 -0.784314 0 | |
| 6.565934 -4.985994 0 | |
| 2.335165 -7.226891 1 | |
| -3.818681 1.792717 0 | |
| 5.412088 -2.801120 0 | |
| 4.862637 -2.016807 0 | |
| -9.038462 4.201681 1 | |
| -9.313187 -0.616246 1 | |
| -0.247253 2.464986 0 | |
| 7.829670 -1.736695 0 | |
| -8.928571 1.064426 1 | |
| 7.719780 -2.857143 0 | |
| 4.752747 -5.042017 0 | |
| -4.917582 -8.907563 1 | |
| -0.851648 -4.425770 1 | |
| -9.258242 -1.344538 1 | |
| -7.994505 -3.081232 1 | |
| -7.774725 0.056022 1 | |
| 6.346154 1.288515 0 | |
| 6.675824 -1.792717 0 | |
| -5.082418 -1.120448 1 | |
| 4.203297 3.361345 0 | |
| 1.071429 2.801120 0 | |
| -8.214286 1.792717 1 | |
| -9.203297 -3.249300 1 | |
| 1.620879 -9.355742 1 | |
| 1.895604 3.473389 0 | |
| 0.631868 -5.994398 1 | |
| 5.027473 -0.280112 0 | |
| -3.708791 -2.745098 1 | |
| -1.620879 -8.123249 1 | |
| 7.170330 -6.442577 0 | |
| -8.598901 -8.011204 1 | |
| 7.829670 0.728291 0 | |
| -1.236264 1.008403 0 | |
| -7.390110 -4.201681 1 | |
| 3.928571 0.784314 0 | |
| -2.060440 -4.089636 1 | |
| 4.368132 7.002801 0 | |
| -0.686813 -6.442577 1 | |
| 5.686813 -8.627451 1 | |
| 4.862637 3.193277 0 | |
| -2.884615 -3.697479 1 | |
| -1.456044 -7.507003 1 | |
| -7.609890 -5.658263 1 | |
| -1.456044 -6.946779 1 | |
| -0.686813 3.305322 0 | |
| -0.302198 -2.577031 0 | |
| 4.093407 -7.787115 1 | |
| -3.653846 -7.170868 1 | |
| -0.247253 4.705882 0 | |
| -1.346154 1.400560 0 | |
| 7.719780 -4.873950 0 | |
| -2.884615 5.994398 0 | |
| 3.324176 5.826331 0 | |
| 7.445055 -0.784314 0 | |
| 0.686813 7.619048 0 | |
| 6.675824 5.882353 0 | |
| 1.126374 -9.355742 1 | |
| 4.478022 -7.787115 1 | |
| -0.906593 3.697479 0 | |
| 4.258242 -0.112045 0 | |
| 6.510989 -0.392157 0 | |
| 6.510989 3.865546 0 | |
| 5.412088 -6.274510 0 | |
| -8.104396 -0.784314 1 | |
| -8.983516 -3.865546 1 | |
| 1.620879 1.512605 0 | |
| 4.862637 1.680672 0 | |
| 6.016484 -4.761905 0 | |
| -5.741758 -1.176471 1 | |
| -6.236264 -0.336134 1 | |
| -7.609890 -0.336134 1 | |
| -1.016484 8.571429 0 | |
| -4.862637 -2.128852 1 | |
| 2.994505 -2.913165 0 | |
| 6.950549 -3.249300 0 | |
| -6.950549 -3.641457 1 | |
| 2.500000 5.882353 0 | |
| -6.401099 -7.170868 1 | |
| -2.280220 2.745098 0 | |
| -2.609890 1.232493 0 | |
| -4.313187 -5.098039 1 | |
| -0.741758 -8.683473 1 | |
| 3.379121 -4.257703 0 | |
| -3.379121 0.448179 0 | |
| -5.961538 -7.394958 1 | |
| 1.510989 -5.490196 1 | |
| 0.247253 1.456583 0 | |
| -7.664835 1.288515 1 | |
| 2.554945 -0.728291 0 | |
| -6.126374 -8.907563 1 | |
| -4.368132 -3.473389 1 | |
| -0.027473 -8.963585 1 | |
| 1.840659 8.123249 0 | |
| 3.049451 4.369748 0 | |
| 2.005495 -5.826331 1 | |
| 6.950549 4.089636 0 | |
| 6.181319 -8.907563 1 | |
| -5.082418 -3.977591 1 | |
| 1.346154 3.137255 0 | |
| 5.631868 -1.176471 0 | |
| 3.159341 -6.666667 1 | |
| -9.093407 1.792717 1 | |
| -5.302198 -1.456583 1 | |
| -7.445055 -8.683473 1 | |
| -4.917582 3.921569 0 | |
| 3.269231 2.857143 0 | |
| -3.269231 -3.081232 1 | |
| 5.412088 5.714286 0 | |
| 2.664835 -2.072829 0 | |
| 1.181319 1.736695 0 | |
| 0.357143 0.224090 0 | |
| -2.225275 -3.305322 1 | |
| 2.280220 2.464986 0 | |
| -2.664835 8.067227 0 | |
| 0.686813 -7.226891 1 | |
| -0.027473 -7.563025 1 | |
| -0.467033 -1.400560 0 | |
| -8.873626 -7.114846 1 | |
| -2.390110 3.417367 0 | |
| 3.983516 -1.400560 0 | |
| 2.335165 -3.193277 0 | |
| -3.928571 -6.162465 1 | |
| -4.038462 -1.904762 1 | |
| 2.060440 0.448179 0 | |
| -2.884615 -7.170868 1 | |
| 3.269231 0.672269 0 | |
| 1.950549 -8.571429 1 | |
| -6.510989 -5.770308 1 | |
| -7.390110 -7.619048 1 | |
| -1.346154 -5.434174 1 | |
| 4.807692 -3.193277 0 | |
| -4.478022 6.330532 0 | |
| -1.181319 6.498599 0 | |
| -8.104396 -8.347339 1 | |
| 5.741758 -9.075630 1 | |
| -3.214286 4.089636 0 | |
| -0.686813 4.761905 0 | |
| 1.016484 -6.442577 1 | |
| 0.247253 -6.722689 1 | |
| -3.434066 -2.352941 1 | |
| 2.500000 -1.736695 0 | |
| -5.796703 -5.770308 1 | |
| -2.994505 7.058824 0 | |
| 0.192308 2.296919 0 | |
| -4.532967 -7.282913 1 | |
| 1.785714 -1.456583 0 | |
| -1.291209 5.602241 0 | |
| -0.686813 -9.355742 1 | |
| 3.104396 -4.873950 0 | |
| 0.906593 -0.504202 0 | |
| -6.565934 -4.649860 1 | |
| -6.565934 -4.649860 1 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| // 必須在演算過程中, 設法限制各變數的上下限! | |
| #include <stdio.h> | |
| #include <stdlib.h> | |
| #include <math.h> | |
| #include <memory.h> | |
| #include <time.h> | |
| // 最大族群數, NP | |
| #define MAXPOP 5000 | |
| // 最大向量維度, D | |
| #define MAXDIM 35 | |
| // MAXIMAPROBLEM =1 最大化 0 最小化 | |
| #define MAXIMAPROBLEM 1 | |
| // 最大化時 PENALITY 必須為負值, 否則為正值 | |
| #define PENALITY -1000 | |
| /* | |
| #define MAXIMAPROBLEM 0 | |
| #define PENALITY 1000 | |
| */ | |
| /*------Constants for rnd_uni()--------------------------------------------*/ | |
| #define IM1 2147483563 | |
| #define IM2 2147483399 | |
| #define AM (1.0/IM1) | |
| #define IMM1 (IM1-1) | |
| #define IA1 40014 | |
| #define IA2 40692 | |
| #define IQ1 53668 | |
| #define IQ2 52774 | |
| #define IR1 12211 | |
| #define IR2 3791 | |
| #define NTAB 32 | |
| #define NDIV (1+IMM1/NTAB) | |
| #define EPS 1.2e-7 | |
| #define RNMX (1.0-EPS) | |
| /*------------------------Globals---------------------------------------*/ | |
| long rnd_uni_init; /* serves as a seed for rnd_uni() */ | |
| double c[MAXPOP][MAXDIM], d[MAXPOP][MAXDIM]; | |
| double (*pold)[MAXPOP][MAXDIM], (*pnew)[MAXPOP][MAXDIM], (*pswap)[MAXPOP][MAXDIM]; | |
| /*---------Function declarations----------------------------------------*/ | |
| void assignd(int D, double a[], double b[]); | |
| double rnd_uni(long *idum); /* uniform pseudo random number generator */ | |
| double extern evaluate(int D, double tmp[], long *nfeval); /* obj. funct. */ | |
| /*---------Function definitions-----------------------------------------*/ | |
| // 指定向量 b 為 a | |
| void assignd(int D, double a[], double b[]) | |
| { | |
| int j; | |
| for (j=0; j<D; j++) | |
| { | |
| a[j] = b[j]; | |
| } | |
| } | |
| // 產生 0 ~ 1 間的亂數 | |
| double rnd_uni(long *idum) | |
| { | |
| long j; | |
| long k; | |
| static long idum2=123456789; | |
| static long iy=0; | |
| static long iv[NTAB]; | |
| double temp; | |
| if (*idum <= 0) | |
| { | |
| if (-(*idum) < 1) *idum=1; | |
| else *idum = -(*idum); | |
| idum2=(*idum); | |
| for (j=NTAB+7;j>=0;j--) | |
| { | |
| k=(*idum)/IQ1; | |
| *idum=IA1*(*idum-k*IQ1)-k*IR1; | |
| if (*idum < 0) *idum += IM1; | |
| if (j < NTAB) iv[j] = *idum; | |
| } | |
| iy=iv[0]; | |
| } | |
| k=(*idum)/IQ1; | |
| *idum=IA1*(*idum-k*IQ1)-k*IR1; | |
| if (*idum < 0) *idum += IM1; | |
| k=idum2/IQ2; | |
| idum2=IA2*(idum2-k*IQ2)-k*IR2; | |
| if (idum2 < 0) idum2 += IM2; | |
| j=iy/NDIV; | |
| iy=iv[j]-idum2; | |
| iv[j] = *idum; | |
| if (iy < 1) iy += IMM1; | |
| if ((temp=AM*iy) > RNMX) return RNMX; | |
| else return temp; | |
| }/*------End of rnd_uni()--------------------------*/ | |
| // 將上下限轉為全域變數 | |
| double inibound_h; /* upper parameter bound */ | |
| double inibound_l; /* lower parameter bound */ | |
| // 與機構合成相關的全域變數 | |
| // 宣告一個座標結構 | |
| struct Coord { | |
| double x; | |
| double y; | |
| // 這裡保留 double z; | |
| }; | |
| main(int argc, char *argv[]) | |
| { | |
| char chr; /* y/n choice variable */ | |
| char *strat[] = /* strategy-indicator */ | |
| { | |
| "", | |
| "DE/best/1/exp", | |
| "DE/rand/1/exp", | |
| "DE/rand-to-best/1/exp", | |
| "DE/best/2/exp", | |
| "DE/rand/2/exp", | |
| "DE/best/1/bin", | |
| "DE/rand/1/bin", | |
| "DE/rand-to-best/1/bin", | |
| "DE/best/2/bin", | |
| "DE/rand/2/bin" | |
| }; | |
| int i, j, L, n; /* counting variables */ | |
| int r1, r2, r3, r4; /* placeholders for random indexes */ | |
| int r5; /* placeholders for random indexes */ | |
| int D; /* Dimension of parameter vector */ | |
| int NP; /* number of population members */ | |
| int imin; /* index to member with lowest energy */ | |
| int refresh; /* refresh rate of screen output */ | |
| int strategy; /* choice parameter for screen output */ | |
| int gen, genmax, seed; | |
| long nfeval; /* number of function evaluations */ | |
| double trial_cost; /* buffer variable */ | |
| // 將上下限轉為全域變數, 可能要根據各變數加以設定 | |
| //double inibound_h; /* upper parameter bound */ | |
| //double inibound_l; /* lower parameter bound */ | |
| double tmp[MAXDIM], best[MAXDIM], bestit[MAXDIM]; /* members */ | |
| double cost[MAXPOP]; /* obj. funct. values */ | |
| double cvar; /* computes the cost variance */ | |
| double cmean; /* mean cost */ | |
| double F,CR; /* control variables of DE */ | |
| double cmin; /* help variables */ | |
| FILE *fpin_ptr; | |
| FILE *fpout_ptr; | |
| // 計算執行過程所需時間起點, 需要導入 time.h | |
| clock_t start = clock(); | |
| /*------Initializations----------------------------*/ | |
| // 將結果寫入 out.dat | |
| fpout_ptr = fopen("out.dat","w"); /* open output file for reading, */ | |
| // 目前已經採用 strategy 3 可以得到最佳結果 | |
| strategy = 3; | |
| genmax = 2000; | |
| refresh = 100; | |
| // 配合機構尺寸合成, 每一個體有 9 個機構尺寸值與 5 個通過點角度值 | |
| D = 2; | |
| NP = 200; | |
| inibound_h = 50.; | |
| inibound_l = 0.; | |
| /*得到最佳解 | |
| F = 0.85; | |
| CR 必須介於 0 to 1. 之間 | |
| CR = 1.; | |
| */ | |
| F = 0.85; | |
| CR = 1.; | |
| seed = 3; | |
| //fclose(fpin_ptr); | |
| /*-----Checking input variables for proper range----------------------------*/ | |
| if (D > MAXDIM) | |
| { | |
| printf("\nError! D=%d > MAXDIM=%d\n",D,MAXDIM); | |
| exit(1); | |
| } | |
| if (D <= 0) | |
| { | |
| printf("\nError! D=%d, should be > 0\n",D); | |
| exit(1); | |
| } | |
| if (NP > MAXPOP) | |
| { | |
| printf("\nError! NP=%d > MAXPOP=%d\n",NP,MAXPOP); | |
| exit(1); | |
| } | |
| if (NP <= 0) | |
| { | |
| printf("\nError! NP=%d, should be > 0\n",NP); | |
| exit(1); | |
| } | |
| if ((CR < 0) || (CR > 1.0)) | |
| { | |
| printf("\nError! CR=%f, should be ex [0,1]\n",CR); | |
| exit(1); | |
| } | |
| if (seed <= 0) | |
| { | |
| printf("\nError! seed=%d, should be > 0\n",seed); | |
| exit(1); | |
| } | |
| if (refresh <= 0) | |
| { | |
| printf("\nError! refresh=%d, should be > 0\n",refresh); | |
| exit(1); | |
| } | |
| if (genmax <= 0) | |
| { | |
| printf("\nError! genmax=%d, should be > 0\n",genmax); | |
| exit(1); | |
| } | |
| if ((strategy < 0) || (strategy > 10)) | |
| { | |
| printf("\nError! strategy=%d, should be ex {1,2,3,4,5,6,7,8,9,10}\n",strategy); | |
| exit(1); | |
| } | |
| if (inibound_h < inibound_l) | |
| { | |
| printf("\nError! inibound_h=%f < inibound_l=%f\n",inibound_h, inibound_l); | |
| exit(1); | |
| } | |
| /*-----Initialize random number generator-----------------------------*/ | |
| rnd_uni_init = -(long)seed; /* initialization of rnd_uni() */ | |
| nfeval = 0; /* reset number of function evaluations */ | |
| /*------Initialization------------------------------------------------*/ | |
| /*------Right now this part is kept fairly simple and just generates--*/ | |
| /*------random numbers in the range [-initfac, +initfac]. You might---*/ | |
| /*------want to extend the init part such that you can initialize-----*/ | |
| /*------each parameter separately.------------------------------------*/ | |
| for (i=0; i<NP; i++) | |
| { | |
| for (j=0; j<D; j++) /* spread initial population members */ | |
| { | |
| c[i][j] = inibound_l + rnd_uni(&rnd_uni_init)*(inibound_h - inibound_l); | |
| } | |
| cost[i] = evaluate(D,c[i],&nfeval); /* obj. funct. value */ | |
| } | |
| cmin = cost[0]; | |
| imin = 0; | |
| for (i=1; i<NP; i++) | |
| { | |
| if(MAXIMAPROBLEM == 1) | |
| { | |
| // 改為最大化 | |
| if (cost[i]>cmin) | |
| { | |
| cmin = cost[i]; | |
| imin = i; | |
| } | |
| } | |
| else | |
| { | |
| // 最小化問題 | |
| if (cost[i]<cmin) | |
| { | |
| cmin = cost[i]; | |
| imin = i; | |
| } | |
| } | |
| } | |
| assignd(D,best,c[imin]); /* save best member ever */ | |
| assignd(D,bestit,c[imin]); /* save best member of generation */ | |
| pold = &c; /* old population (generation G) */ | |
| pnew = &d; /* new population (generation G+1) */ | |
| /*=======================================================================*/ | |
| /*=========Iteration loop================================================*/ | |
| /*=======================================================================*/ | |
| gen = 0; /* generation counter reset */ | |
| while ((gen < genmax) /*&& (kbhit() == 0)*/) /* remove comments if conio.h */ | |
| { /* is accepted by compiler */ | |
| gen++; | |
| imin = 0; | |
| for (i=0; i<NP; i++) /* Start of loop through ensemble */ | |
| { | |
| do /* Pick a random population member */ | |
| { /* Endless loop for NP < 2 !!! */ | |
| r1 = (int)(rnd_uni(&rnd_uni_init)*NP); | |
| }while(r1==i); | |
| do /* Pick a random population member */ | |
| { /* Endless loop for NP < 3 !!! */ | |
| r2 = (int)(rnd_uni(&rnd_uni_init)*NP); | |
| }while((r2==i) || (r2==r1)); | |
| do /* Pick a random population member */ | |
| { /* Endless loop for NP < 4 !!! */ | |
| r3 = (int)(rnd_uni(&rnd_uni_init)*NP); | |
| }while((r3==i) || (r3==r1) || (r3==r2)); | |
| do /* Pick a random population member */ | |
| { /* Endless loop for NP < 5 !!! */ | |
| r4 = (int)(rnd_uni(&rnd_uni_init)*NP); | |
| }while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3)); | |
| do /* Pick a random population member */ | |
| { /* Endless loop for NP < 6 !!! */ | |
| r5 = (int)(rnd_uni(&rnd_uni_init)*NP); | |
| }while((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4)); | |
| /*=======EXPONENTIAL CROSSOVER============================================================*/ | |
| /*-------DE/best/1/exp--------------------------------------------------------------------*/ | |
| /*-------Our oldest strategy but still not bad. However, we have found several------------*/ | |
| /*-------optimization problems where misconvergence occurs.-------------------------------*/ | |
| if (strategy == 1) /* strategy DE0 (not in our paper) */ | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| L = 0; | |
| do | |
| { | |
| tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]); | |
| n = (n+1)%D; | |
| L++; | |
| }while((rnd_uni(&rnd_uni_init) < CR) && (L < D)); | |
| } | |
| /*-------DE/rand/1/exp-------------------------------------------------------------------*/ | |
| /*-------This is one of my favourite strategies. It works especially well when the-------*/ | |
| /*-------"bestit[]"-schemes experience misconvergence. Try e.g. F=0.7 and CR=0.5---------*/ | |
| /*-------as a first guess.---------------------------------------------------------------*/ | |
| else if (strategy == 2) /* strategy DE1 in the techreport */ | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| L = 0; | |
| do | |
| { | |
| tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]); | |
| n = (n+1)%D; | |
| L++; | |
| }while((rnd_uni(&rnd_uni_init) < CR) && (L < D)); | |
| } | |
| /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/ | |
| /*-------This strategy seems to be one of the best strategies. Try F=0.85 and CR=1.------*/ | |
| /*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/ | |
| /*-------should play around with all three control variables.----------------------------*/ | |
| else if (strategy == 3) /* similiar to DE2 but generally better */ | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| L = 0; | |
| do | |
| { | |
| tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]); | |
| n = (n+1)%D; | |
| L++; | |
| }while((rnd_uni(&rnd_uni_init) < CR) && (L < D)); | |
| } | |
| /*-------DE/best/2/exp is another powerful strategy worth trying--------------------------*/ | |
| else if (strategy == 4) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| L = 0; | |
| do | |
| { | |
| tmp[n] = bestit[n] + | |
| ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F; | |
| n = (n+1)%D; | |
| L++; | |
| }while((rnd_uni(&rnd_uni_init) < CR) && (L < D)); | |
| } | |
| /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/ | |
| else if (strategy == 5) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| L = 0; | |
| do | |
| { | |
| tmp[n] = (*pold)[r5][n] + | |
| ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F; | |
| n = (n+1)%D; | |
| L++; | |
| }while((rnd_uni(&rnd_uni_init) < CR) && (L < D)); | |
| } | |
| /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/ | |
| /*-------DE/best/1/bin--------------------------------------------------------------------*/ | |
| else if (strategy == 6) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| for (L=0; L<D; L++) /* perform D binomial trials */ | |
| { | |
| if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */ | |
| { | |
| tmp[n] = bestit[n] + F*((*pold)[r2][n]-(*pold)[r3][n]); | |
| } | |
| n = (n+1)%D; | |
| } | |
| } | |
| /*-------DE/rand/1/bin-------------------------------------------------------------------*/ | |
| else if (strategy == 7) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| for (L=0; L<D; L++) /* perform D binomial trials */ | |
| { | |
| if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */ | |
| { | |
| tmp[n] = (*pold)[r1][n] + F*((*pold)[r2][n]-(*pold)[r3][n]); | |
| } | |
| n = (n+1)%D; | |
| } | |
| } | |
| /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/ | |
| else if (strategy == 8) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| for (L=0; L<D; L++) /* perform D binomial trials */ | |
| { | |
| if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */ | |
| { | |
| tmp[n] = tmp[n] + F*(bestit[n] - tmp[n]) + F*((*pold)[r1][n]-(*pold)[r2][n]); | |
| } | |
| n = (n+1)%D; | |
| } | |
| } | |
| /*-------DE/best/2/bin--------------------------------------------------------------------*/ | |
| else if (strategy == 9) | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| for (L=0; L<D; L++) /* perform D binomial trials */ | |
| { | |
| if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */ | |
| { | |
| tmp[n] = bestit[n] + | |
| ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F; | |
| } | |
| n = (n+1)%D; | |
| } | |
| } | |
| /*-------DE/rand/2/bin--------------------------------------------------------------------*/ | |
| else | |
| { | |
| assignd(D,tmp,(*pold)[i]); | |
| n = (int)(rnd_uni(&rnd_uni_init)*D); | |
| for (L=0; L<D; L++) /* perform D binomial trials */ | |
| { | |
| if ((rnd_uni(&rnd_uni_init) < CR) || L == (D-1)) /* change at least one parameter */ | |
| { | |
| tmp[n] = (*pold)[r5][n] + | |
| ((*pold)[r1][n]+(*pold)[r2][n]-(*pold)[r3][n]-(*pold)[r4][n])*F; | |
| } | |
| n = (n+1)%D; | |
| } | |
| } | |
| /*=======Trial mutation now in tmp[]. Test how good this choice really was.==================*/ | |
| trial_cost = evaluate(D,tmp,&nfeval); /* Evaluate new vector in tmp[] */ | |
| if(MAXIMAPROBLEM == 1) | |
| { | |
| // 改為最大化 | |
| if (trial_cost >= cost[i]) /* improved objective function value ? */ | |
| { | |
| cost[i]=trial_cost; | |
| assignd(D,(*pnew)[i],tmp); | |
| if (trial_cost>cmin) /* Was this a new minimum? */ | |
| { /* if so...*/ | |
| cmin=trial_cost; /* reset cmin to new low...*/ | |
| imin=i; | |
| assignd(D,best,tmp); | |
| } | |
| } | |
| else | |
| { | |
| assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */ | |
| } | |
| } | |
| else | |
| { | |
| // 最小化問題 | |
| if (trial_cost <= cost[i]) /* improved objective function value ? */ | |
| { | |
| cost[i]=trial_cost; | |
| assignd(D,(*pnew)[i],tmp); | |
| if (trial_cost<cmin) /* Was this a new minimum? */ | |
| { /* if so...*/ | |
| cmin=trial_cost; /* reset cmin to new low...*/ | |
| imin=i; | |
| assignd(D,best,tmp); | |
| } | |
| } | |
| else | |
| { | |
| assignd(D,(*pnew)[i],(*pold)[i]); /* replace target with old value */ | |
| } | |
| } | |
| } /* End mutation loop through pop. */ | |
| assignd(D,bestit,best); /* Save best population member of current iteration */ | |
| /* swap population arrays. New generation becomes old one */ | |
| pswap = pold; | |
| pold = pnew; | |
| pnew = pswap; | |
| /*----Compute the energy variance (just for monitoring purposes)-----------*/ | |
| cmean = 0.; /* compute the mean value first */ | |
| for (j=0; j<NP; j++) | |
| { | |
| cmean += cost[j]; | |
| } | |
| cmean = cmean/NP; | |
| cvar = 0.; /* now the variance */ | |
| for (j=0; j<NP; j++) | |
| { | |
| cvar += (cost[j] - cmean)*(cost[j] - cmean); | |
| } | |
| cvar = cvar/(NP-1); | |
| /*----Output part----------------------------------------------------------*/ | |
| if (gen%refresh==1) /* display after every refresh generations */ | |
| { /* ABORT works only if conio.h is accepted by your compiler */ | |
| printf("\n\n PRESS ANY KEY TO ABORT"); | |
| printf("\n\n\n Best-so-far cost funct. value=%-15.10g\n",cmin); | |
| for (j=0;j<D;j++) | |
| { | |
| printf("\n best[%d]=%-15.10g",j,best[j]); | |
| } | |
| printf("\n\n Generation=%d NFEs=%ld Strategy: %s ",gen,nfeval,strat[strategy]); | |
| printf("\n NP=%d F=%-4.2g CR=%-4.2g cost-variance=%-10.5g\n", | |
| NP,F,CR,cvar); | |
| } | |
| fprintf(fpout_ptr,"%ld %-15.10g\n",nfeval,cmin); | |
| } | |
| /*=======================================================================*/ | |
| /*=========End of iteration loop=========================================*/ | |
| /*=======================================================================*/ | |
| /*-------Final output in file-------------------------------------------*/ | |
| fprintf(fpout_ptr,"\n\n\n Best-so-far obj. funct. value = %-15.10g\n",cmin); | |
| for (j=0;j<D;j++) | |
| { | |
| fprintf(fpout_ptr,"\n best[%d]=%-15.10g",j,best[j]); | |
| } | |
| fprintf(fpout_ptr,"\n\n Generation=%d NFEs=%ld Strategy: %s ",gen,nfeval,strat[strategy]); | |
| fprintf(fpout_ptr,"\n NP=%d F=%-4.2g CR=%-4.2g cost-variance=%-10.5g\n", | |
| NP,F,CR,cvar); | |
| fclose(fpout_ptr); | |
| /* Code you want timed here */ | |
| printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC); | |
| return(0); | |
| } | |
| /*-----------End of main()------------------------------------------*/ | |
| // 適應函式 fittness function (cost function) | |
| double evaluate(int D, double tmp[], long *nfeval) | |
| { | |
| double result=0, surface = 192.0, z, volume, penality; | |
| (*nfeval)++; | |
| z = (surface-tmp[0]*tmp[1])/(2.0*(tmp[0]+tmp[1])); | |
| volume = tmp[0]*tmp[1]*z; | |
| if(volume <= 0){ | |
| return PENALITY; | |
| } | |
| // 只限制長度與寬度必須大於 0 | |
| if(tmp[0] <= inibound_l){ | |
| return PENALITY; | |
| } | |
| if(tmp[1] <= inibound_l){ | |
| return PENALITY; | |
| } | |
| /* | |
| if((tmp[0] <= inibound_l)|| (tmp[0] >inibound_h)){ | |
| return PENALITY; | |
| } | |
| if((tmp[1] <= inibound_l) || (tmp[1] >inibound_h)){ | |
| return PENALITY; | |
| } | |
| */ | |
| // volume must >0 and max volume | |
| // 目前為最小化問題 | |
| //return 1+1/(volume*volume); | |
| return volume; | |
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import random | |
| class Chromosome(object): | |
| """ | |
| just copy the idea of genetic algorithm, pretty similar.. | |
| """ | |
| def __init__(self, n): | |
| """ | |
| int n, dimension of question | |
| """ | |
| # dimension | |
| self.n = n | |
| # the gene | |
| self.v = [0] * n | |
| # the fitness value | |
| self.f = 0 | |
| def assign(self, obj): | |
| """ | |
| Chromosome obj | |
| copy all attribute from obj to itself | |
| """ | |
| self.n = obj.n | |
| self.v = obj.v[:] | |
| self.f = obj.f | |
| class DiffertialEvolution(object): | |
| def __init__(self,Func, pType, strategy, D, NP, F, CR, lower, upper, maxGen, report): | |
| # if pType = 1 it is maximization otherwise is minimization problem | |
| self.pType = pType | |
| # strategy 1~10, choice what strategy to generate new member in temporary | |
| self.strategy = strategy | |
| # dimesion of quesiton | |
| self.D = D | |
| # population size | |
| # To start off NP = 10*D is a reasonable choice. Increase NP if misconvergence | |
| self.NP = NP | |
| # weight factor | |
| # F is usually between 0.5 and 1 (in rare cases > 1) | |
| self.F = F | |
| # crossover possible | |
| # CR in [0,1] | |
| self.CR = CR | |
| # lower bound array | |
| self.lb = lower[:] | |
| # upper bound array | |
| self.ub = upper[:] | |
| # maximum generation | |
| self.maxGen = maxGen | |
| # how many generation report once | |
| self.rpt = report | |
| # object function, or enviorment | |
| self.f = Func | |
| # check parameter is set properly | |
| self.checkParameter() | |
| # generation pool, depend on population size | |
| self.pop = [Chromosome(D) for i in range(NP)] | |
| # last generation best member | |
| self.lastgenbest = Chromosome(D) | |
| # current best member | |
| self.currentbest = Chromosome(D) | |
| # the generation count | |
| self.gen = 0 | |
| # the vector | |
| self.r1 = 0 | |
| self.r2 = 0 | |
| self.r3 = 0 | |
| self.r4 = 0 | |
| self.r5 = 0 | |
| def checkParameter(self): | |
| """ | |
| check parameter is set properly | |
| """ | |
| if (type(self.D) is not int) and self.D <= 0: | |
| raise Exception('D shoud be integer and larger than 0') | |
| if (type(self.NP) is not int) and self.NP <= 0: | |
| raise Exception('NP shoud be integer and larger than 0') | |
| if self.CR < 0 or self.CR > 1: | |
| raise Exception('CR should be [0,1]') | |
| if self.maxGen <= 0: | |
| raise Exception('generation should larger than 0') | |
| if self.rpt <= 0 or self.rpt > self.maxGen: | |
| raise Exception('report should be larger than 0 and less than max genration') | |
| if self.strategy < 1 or self.strategy > 10: | |
| raise Exception('strategy should be [1,10]') | |
| for lower, upper in zip(self.lb, self.ub): | |
| if lower > upper: | |
| raise Exception('upper bound should be larger than lower bound') | |
| def init(self): | |
| """ | |
| init population | |
| """ | |
| for i in range(self.NP): | |
| for j in range(self.D): | |
| self.pop[i].v[j] = self.lb[j] + random.random()*(self.ub[j] - self.lb[j]) | |
| self.pop[i].f = self.evalute(self.pop[i]) | |
| def evalute(self, p): | |
| """ | |
| evalute the member in enviorment | |
| """ | |
| return self.f(p.v) | |
| def findBest(self): | |
| """ | |
| find member that have minimum fitness value from pool | |
| """ | |
| if self.pType == 1: | |
| return max(self.pop, key=lambda chrom:chrom.f) | |
| else: | |
| return min(self.pop, key=lambda chrom:chrom.f) | |
| def generateRandomVector(self, i): | |
| """ | |
| generate new vector | |
| """ | |
| while True: | |
| self.r1 = int(random.random() * self.NP) | |
| if not (self.r1 == i): | |
| break | |
| while True: | |
| self.r2 = int(random.random() * self.NP) | |
| if not ((self.r2 == i) or (self.r2 == self.r1)): | |
| break | |
| while True: | |
| self.r3 = int(random.random() * self.NP) | |
| if not ((self.r3 == i) or (self.r3 == self.r1) or (self.r3 == self.r2)): | |
| break | |
| while True: | |
| self.r4 = int(random.random() * self.NP) | |
| if not ((self.r4 == i) or (self.r4 == self.r1) or (self.r4 == self.r2) or (self.r4 == self.r3)): | |
| break | |
| while True: | |
| self.r5 = int(random.random() * self.NP) | |
| if not ((self.r5 == i) or (self.r5 == self.r1) or (self.r5 == self.r2) or (self.r5 == self.r3) or (self.r5 == self.r4)): | |
| break | |
| def recombination(self, i): | |
| """ | |
| use new vector, recombination the new one member to tmp | |
| """ | |
| tmp = Chromosome(self.D) | |
| if self.strategy == 1: | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| L = 0 | |
| while True: | |
| tmp.v[n] = self.lastgenbest.v[n] + self.F*(self.pop[self.r2].v[n] - self.pop[self.r3].v[n]) | |
| n = (n + 1) % self.D | |
| L += 1 | |
| if not ((random.random() < self.CR) and (L < self.D)): | |
| break | |
| elif self.strategy == 2: | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| L = 0 | |
| while True: | |
| tmp.v[n] = self.pop[self.r1].v[n] + self.F*(self.pop[self.r2].v[n] - self.pop[self.r3].v[n]) | |
| n = (n + 1) % self.D | |
| L += 1 | |
| if not ((random.random() < self.CR) and (L < self.D)): | |
| break | |
| elif (self.strategy == 3): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| L = 0 | |
| while True: | |
| tmp.v[n] = tmp.v[n] + self.F*(self.lastgenbest.v[n] - tmp.v[n]) + self.F*(self.pop[self.r1].v[n] - self.pop[self.r2].v[n]) | |
| n = (n + 1) % self.D | |
| L += 1 | |
| if not ((random.random() < self.CR) and (L < self.D)): | |
| break | |
| elif (self.strategy == 4): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| L = 0 | |
| while True: | |
| tmp.v[n] = self.lastgenbest.v[n] + (self.pop[self.r1].v[n] + self.pop[self.r2].v[n] - self.pop[self.r3].v[n] - self.pop[self.r4].v[n]) * self.F | |
| n = (n + 1) % self.D | |
| L += 1 | |
| if not ((random.random() < self.CR) and (L < self.D)): | |
| break | |
| elif (self.strategy == 5): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| L = 0 | |
| while True: | |
| tmp.v[n] = self.pop[self.r5].v[n] + (self.pop[self.r1].v[n] + self.pop[self.r2].v[n] - self.pop[self.r3].v[n] - self.pop[self.r4].v[n]) * self.F | |
| n = (n + 1) % self.D | |
| L += 1 | |
| if not ((random.random() < self.CR) and (L < self.D)): | |
| break | |
| elif (self.strategy == 6): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| for L in range(self.D): | |
| if ((random.random() < self.CR) or L == (self.D - 1)): | |
| tmp.v[n] = self.lastgenbest.v[n] + self.F*(self.pop[self.r2].v[n] - self.pop[self.r3].v[n]) | |
| n = (n + 1) % self.D | |
| elif (self.strategy == 7): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| for L in range(self.D): | |
| if ((random.random() < self.CR) or L == (self.D - 1)): | |
| tmp.v[n] = self.pop[self.r1].v[n] + self.F*(self.pop[self.r2].v[n] - self.pop[self.r3].v[n]) | |
| n = (n + 1) % self.D | |
| elif (self.strategy == 8): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| for L in range(self.D): | |
| if ((random.random() < self.CR) or L == (self.D - 1)): | |
| tmp.v[n] = tmp.v[n] + self.F*(self.lastgenbest.v[n] - tmp.v[n]) + self.F*(self.pop[self.r1].v[n] - self.pop[self.r2].v[n]) | |
| n = (n + 1) % self.D | |
| elif (self.strategy == 9): | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| for L in range(self.D): | |
| if ((random.random() < self.CR) or L == (self.D - 1)): | |
| tmp.v[n] = self.lastgenbest.v[n] + (self.pop[self.r1].v[n] + self.pop[self.r2].v[n] - self.pop[self.r3].v[n] - self.pop[self.r4].v[n]) * self.F | |
| n = (n + 1) % self.D | |
| else: | |
| tmp.assign(self.pop[i]) | |
| n = int(random.random() * self.D) | |
| for L in range(self.D): | |
| if ((random.random() < self.CR) or L == (self.D - 1)): | |
| tmp.v[n] = self.pop[self.r5].v[n] + (self.pop[self.r1].v[n] + self.pop[self.r2].v[n] - self.pop[self.r3].v[n] - self.pop[self.r4].v[n]) * self.F | |
| n = (n + 1) % self.D | |
| return tmp | |
| def report(self): | |
| """ | |
| report current generation status | |
| """ | |
| if self.gen == 0: | |
| print("DiffertialEvolution results - init pop") | |
| elif self.gen == self.maxGen: | |
| print("Final DiffertialEvolution results at", self.gen, "generations") | |
| else: | |
| print("DiffertialEvolution results after", self.gen, "generations") | |
| print("Function : %.6f" % (self.currentbest.f)) | |
| for i, v in enumerate(self.currentbest.v, start=1): | |
| print("Var", i, ":", v) | |
| def overbound(self, member): | |
| """ | |
| check the member's chromosome that is out of bound? | |
| """ | |
| for i in range(self.D): | |
| if member.v[i] > self.ub[i] or member.v[i] < self.lb[i]: | |
| return True | |
| return False | |
| def run(self): | |
| """ | |
| run the algorithm... | |
| """ | |
| # initial step | |
| # generation 0 | |
| self.gen = 0 | |
| # init the member's chromsome | |
| self.init() | |
| # find the best one(smallest fitness value) | |
| tmp = self.findBest() | |
| # copy to lastgenbest | |
| self.lastgenbest.assign(tmp) | |
| # copy to currentbest | |
| self.currentbest.assign(tmp) | |
| # report status | |
| self.report() | |
| # end initial step | |
| # the evolution journey is beggin... | |
| for self.gen in range(1, self.maxGen + 1): | |
| for i in range(self.NP): | |
| # generate new vector | |
| self.generateRandomVector(i) | |
| # use the vector recombine the member to temporary | |
| tmp = self.recombination(i) | |
| # check the one is out of bound? | |
| if self.overbound(tmp): | |
| # if it is, then ignore | |
| continue | |
| # is not out of bound, that mean it's quilify of enviorment | |
| # then evalute the one | |
| tmp.f = self.evalute(tmp) | |
| # if temporary one is better than origin(fitness value is larger or smaller) | |
| # pType is 1, the problem is maximization type | |
| if self.pType == 1: | |
| if tmp.f >= self.pop[i].f: | |
| # copy the temporary one to origin member | |
| self.pop[i].assign(tmp) | |
| # check the temporary one is better than the currentbest | |
| if tmp.f > self.currentbest.f: | |
| # copy the temporary one to currentbest | |
| self.currentbest.assign(tmp) | |
| else: | |
| if tmp.f <= self.pop[i].f: | |
| # copy the temporary one to origin member | |
| self.pop[i].assign(tmp) | |
| # check the temporary one is better than the currentbest | |
| if tmp.f < self.currentbest.f: | |
| # copy the temporary one to currentbest | |
| self.currentbest.assign(tmp) | |
| # copy the currentbest to lastgenbest | |
| self.lastgenbest.assign(self.currentbest) | |
| # if report generation is set, report | |
| if self.rpt != 0: | |
| if self.gen % self.rpt == 0: | |
| self.report() | |
| # the evolution journey is done, report the final status | |
| self.report() | |
| #fittness function (cost function) | |
| def evaluate(designVariablel): | |
| surface = 80.0 | |
| # if pType is 1, the penality is negative (maximization problem) | |
| # if pType is 0, the penality is positive (minimization problem) | |
| penality = -1000 | |
| z = (surface-designVariablel[0]*designVariablel[1])/(2.0*(designVariablel[0]\ | |
| +designVariablel[1])) | |
| volume = designVariablel[0]*designVariablel[1]*z | |
| if(volume <= 0): | |
| return penality | |
| # box length and width need to be larger than 0 | |
| if(designVariablel[0] <= 0): | |
| return penality | |
| if(designVariablel[1] <= 0): | |
| return penality | |
| return volume | |
| #volume = DiffertialEvolution((self,Func, max, strategy, D, NP, F, CR, lower, upper, maxGen, report) | |
| volume = DiffertialEvolution(evaluate, 1, 3, 2, 100, 0.6, 0.85, [0, 0], [50, 50], 100, 10) | |
| volume.run() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment