# include <stdlib.h>
# include <stdio.h>
# include <string.h>
# include <math.h>
# include <time.h>
# include "fem1d_lagrange.h"
/******************************************************************************/
void fem1d_lagrange_stiffness ( int x_num, double x[], int q_num,
double f ( double x ), double a[], double m[], double b[] )
/******************************************************************************/
/*
Purpose:
FEM1D_LAGRANGE_STIFFNESS evaluates the Lagrange polynomial stiffness matrix.
Discussion:
The finite element method is to be applied over a given interval that
has been meshed with X_NUM points X.
The finite element basis functions are to be the X_NUM Lagrange
basis polynomials L(i)(X), such that
L(i)(X(j)) = delta(i,j).
The following items are computed:
* A, the stiffness matrix, with A(I,J) = integral L'(i)(x) L'(j)(x)
* M, the mass matrix, with M(I,J) = integral L(i)(x) L(j)(x)
* B, the load matrix, with B(I) = integral L(i)(x) F(x)
The integrals are approximated by quadrature.
Boundary conditions are not handled by this routine.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 November 2014
Author:
John Burkardt
Parameters:
Input, int X_NUM, the number of nodes.
Input, double X[X_NUM], the coordinates of the nodes.
Input, int Q_NUM, the number of quadrature points to use.
Input, double F ( double X ), the right hand side function.
Output, double A[X_NUM*X_NUM], the stiffness matrix.
Output, double M[X_NUM*X_NUM], the mass matrix.
Output, double B[X_NUM], the right hand side vector.
*/
{
double *l;
double li;
double lj;
double *lp;
double lpi;
double lpj;
int q_i;
double *q_w;
double *q_x;
int x_i;
int x_j;
/*
Get the quadrature rule for [-1,+1].
*/
q_x = ( double * ) malloc ( q_num * sizeof ( double ) );
q_w = ( double * ) malloc ( q_num * sizeof ( double ) );
legendre_set ( q_num, q_x, q_w );
/*
Adjust the quadrature rule to the interval [ x(1), x(x_num) }
*/
for ( q_i = 0; q_i < q_num; q_i++ )
{
q_x[q_i] = ( ( 1.0 - q_x[q_i] ) * x[0]
+ ( 1.0 + q_x[q_i] ) * x[x_num-1] )
/ 2.0;
q_w[q_i] = q_w[q_i] * ( x[x_num-1] - x[0] ) / 2.0;
}
/*
Evaluate all the Lagrange basis polynomials at all the quadrature points.
*/
l = lagrange_value ( x_num, x, q_num, q_x );
/*
Evaluate all the Lagrange basis polynomial derivatives at all
the quadrature points.
*/
lp = lagrange_derivative ( x_num, x, q_num, q_x );
/*
Assemble the matrix and right hand side.
*/
for ( x_j = 0; x_j < x_num; x_j++ )
{
for ( x_i = 0; x_i < x_num; x_i++ )
{
a[x_i+x_j*x_num] = 0.0;
m[x_i+x_j*x_num] = 0.0;
}
b[x_j] = 0.0;
}
for ( x_i = 0; x_i < x_num; x_i++ )
{
for ( q_i = 0; q_i < q_num; q_i++ )
{
li = l[q_i+x_i*q_num];
lpi = lp[q_i+x_i*q_num];
for ( x_j = 0; x_j < x_num; x_j++ )
{
lj = l[q_i+x_j*q_num];
lpj = lp[q_i+x_j*q_num];
a[x_i+x_j*x_num] = a[x_i+x_j*x_num] + q_w[q_i] * lpi * lpj;
m[x_i+x_j*x_num] = m[x_i+x_j*x_num] + q_w[q_i] * li * lj;
}
b[x_i] = b[x_i] + q_w[q_i] * li * f ( q_x[q_i] );
}
}
free ( l );
free ( lp );
free ( q_w );
free ( q_x );
return;
}
/******************************************************************************/
double *lagrange_derivative ( int nd, double xd[], int ni, double xi[] )
/******************************************************************************/
/*
Purpose:
LAGRANGE_DERIVATIVE evaluates the Lagrange basis derivative.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
19 November 2014
Author:
John Burkardt
Parameters:
Input, int ND, the number of data points.
ND must be at least 1.
Input, double XD[ND], the data points.
Input, int NI, the number of interpolation points.
Input, double XI[NI], the interpolation points.
Output, double LAGRANGE_DERIVATIVE[NI*ND], the values
of the Lagrange basis derivatives at the interpolation points.
*/
{
int i;
int j;
int j1;
int j2;
double *lpi;
double p;
lpi = ( double * ) malloc ( ni * nd * sizeof ( double ) );
for ( j = 0; j < nd; j++ )
{
for ( i = 0; i < ni; i++ )
{
lpi[i+j*ni] = 0.0;
}
}
for ( i = 0; i < ni; i++ )
{
for ( j = 0; j < nd; j++ )
{
for ( j1 = 0; j1 < nd; j1++ )
{
if ( j1 != j )
{
p = 1.0;
for ( j2 = 0; j2 < nd; j2++ )
{
if ( j2 == j1 )
{
p = p / ( xd[j] - xd[j2] );
}
else if ( j2 != j )
{
p = p * ( xi[i] - xd[j2] ) / ( xd[j] - xd[j2] );
}
}
lpi[i+j*ni] = lpi[i+j*ni] + p;
}
}
}
}
return lpi;
}
/******************************************************************************/
double *lagrange_value ( int nd, double xd[], int ni, double xi[] )
/******************************************************************************/
/*
Purpose:
LAGRANGE_VALUE evaluates the Lagrange basis polynomials.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
13 September 2012
Author:
John Burkardt
Parameters:
Input, int ND, the number of data points.
ND must be at least 1.
Input, double XD[ND], the data points.
Input, int NI, the number of interpolation points.
Input, double XI[NI], the interpolation points.
Output, double LAGRANGE_BASIS[NI*ND], the values
of the Lagrange basis polynomials at the interpolation points.
*/
{
int i;
int j;
int k;
double *li;
/*
Evaluate the polynomial.
*/
li = ( double * ) malloc ( ni * nd * sizeof ( double ) );
for ( j = 0; j < nd; j++ )
{
for ( i = 0; i < ni; i++ )
{
li[i+j*ni] = 1.0;
}
}
for ( i = 0; i < nd; i++ )
{
for ( j = 0; j < nd; j++ )
{
if ( j != i )
{
for ( k = 0; k < ni; k++ )
{
li[k+i*ni] = li[k+i*ni] * ( xi[k] - xd[j] ) / ( xd[i] - xd[j] );
}
}
}
}
return li;
}
/******************************************************************************/
void legendre_set ( int n, double x[], double w[] )
/******************************************************************************/
/*
Purpose:
LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature.
Discussion:
The integral:
Integral ( -1 <= X <= 1 ) F(X) dX
Quadrature rule:
Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
The quadrature rule is exact for polynomials through degree 2*N-1.
The abscissas are the zeroes of the Legendre polynomial P(ORDER)(X).
Mathematica can compute the abscissas and weights of a Gauss-Legendre
rule of order N for the interval [A,B] with P digits of precision
by the commands:
Needs["NumericalDifferentialEquationAnalysis`"]
GaussianQuadratureWeights [n, a, b, p ]
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
20 April 2010
Author:
John Burkardt
Reference:
Milton Abramowitz, Irene Stegun,
Handbook of Mathematical Functions,
National Bureau of Standards, 1964,
ISBN: 0-486-61272-4,
LC: QA47.A34.
Vladimir Krylov,
Approximate Calculation of Integrals,
Dover, 2006,
ISBN: 0486445798,
LC: QA311.K713.
Arthur Stroud, Don Secrest,
Gaussian Quadrature Formulas,
Prentice Hall, 1966,
LC: QA299.4G3S7.
Stephen Wolfram,
The Mathematica Book,
Fourth Edition,
Cambridge University Press, 1999,
ISBN: 0-521-64314-7,
LC: QA76.95.W65.
Daniel Zwillinger, editor,
CRC Standard Mathematical Tables and Formulae,
30th Edition,
CRC Pr
C 代码 设置与有限元关联的矩阵和向量 边界值问题 (BVP) 的方法 (FEM).rar
版权申诉
179 浏览量
2022-11-13
21:16:25
上传
评论
收藏 8KB RAR 举报
卷积神经网络
- 粉丝: 338
- 资源: 8460