# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
int main ( int argc, char *argv[] );
void assemble_newton ( double adiag[], double aleft[], double arite[],
double f[], double fold[], double h[], int indx[], int n, int nl,
int node[], int nquad, int nu, int problem, double ul, double ur,
double xn[], double xquad[] );
void assemble_picard ( double adiag[], double aleft[], double arite[],
double f[], double fold[], double h[], int indx[], int n, int nl,
int node[], int nquad, int nu, int problem, double ul, double ur,
double xn[], double xquad[] );
void compare ( double f[], int indx[], int n, int nl, int node[], int nprint,
int nu, int problem, double ul, double ur, double xl, double xn[],
double xr );
double ff ( double x, int problem );
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 *imax, int *nprint, int *nquad, int problem,
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, int problem );
void prsys ( double adiag[], double aleft[], double arite[], double f[],
int nu );
double qq ( double x, int problem );
void solve ( double adiag[], double aleft[], double arite[], double f[],
int nu );
void timestamp ( );
double u_exact ( double x, int problem );
/******************************************************************************/
int main ( int argc, char *argv[] )
/******************************************************************************/
/*
Purpose:
MAIN is the main program for FEM1D_NONLINEAR.
Discussion:
FEM1D_NONLINLEAR solves a nonlinear one dimensional boundary value problem.
The differential equation has the form:
-d/dx (p(x) du/dx) + q(x)*u + u*u' = f(x)
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.
Sample problem #1:
u(x) = x,
p(x) = 1,
q(x) = 0,
f(x) = x,
u(0) = 0,
u'(1) = 1.
The code should solve this problem exactly.
Sample problem #2:
u(x) = 2*(1-cos(0.5*pi*x))/pi,
p(x) = 1,
q(x) = 0,
f(x) = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
u(0) = 0,
u'(1) = 1.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
23 June 2019
Author:
C version by John Burkardt
Parameters:
double ADIAG(NU).
ADIAG(I) is the "diagonal" coefficient of the I-th
equation in the linear system. That is, ADIAG(I) is
the coefficient of the I-th unknown in the I-th equation.
double ALEFT(NU).
ALEFT(I) is the "left hand" coefficient of the I-th
equation in the linear system. 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 FOLD(NU).
FOLD contains the value of F from the previous iteration,
and is used in ASSEMBLE to add correction terms to the
matrix and right hand side.
double H(N), the length of 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 IMAX.
The number of Newton iterations to carry out.
int INDX(0:N).
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 N, which will be unknown N or N+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,N).
For each subinterval I:
NODE(1,I) is the number of the left node, and
NODE(2,I) is the number of the right node.
int NPRINT.
The number of points at which the computed solution
should be printed out when compared to the exact solution.
int NQUAD.
The number of quadrature points used in a subinterval.
This code uses NQUAD = 1.
int N.
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 N-1,
N, or N+1 unknown values, which are the coefficients
of basis functions.
int PROBLEM, indicates which problem to be solved.
* 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
* 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
u(0) = 0, u'(1)=1.
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:N).
XN(I) is the location of the I-th node. XN(0) is XL,
and XN(N) is XR.
double XQUAD(N)
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 N 10
# define NL 2
double adiag[N+1];
double aleft[N+1];
double arite[N+1];
double f[N+1];
double fold[N+1];
double h[N];
int i;
int ibc;
int imax;
int indx[N+1];
int j;
int node[NL*N];
int nprint;
int nquad;
int nu;
int problem;
double ul;
double ur;
double xl;
double xn[N+1];
double xquad[N];
double xr;
timestamp ( );
printf ( "\n" );
printf ( "FEM1D_NONLINEAR\n" );
printf ( " C version\n" );
printf ( "\n" );
printf ( " Solve a nonlinear boundary value problem:\n" );
printf ( "\n" );
printf ( " -d/dx (p(x) du/dx) + q(x)*u + u*u' = f(x)\n" );
printf ( "\n" );
printf ( " on an interval [xl,xr], with the values of\n" );
printf ( " u or u' specified at xl and xr.\n" );
printf ( "\n" );
printf ( " The interval [XL,XR] is broken into N = %d subintervals\n", N