Loading

This lesson concerns a simple one-dimensional two-point boundary value problem. The whole program can be found in Example 1 in the examples directory. We will examine it here below line by line.

- We start by including the header file OFELI.h that itself includes all kernel
class definitions.

#include "OFELI.h"

- The OFELI library is embedded in namespace called
OFELI.

using namespace OFELI;

- Our program has arguments that will described later.

int main(int argc, char *argv[]) {

- L is the length
of the interval in which the problem is defined. Here we have fixed their respective values
at 0 and 1.
N is the number of finite elements. Its default value is 10

double L=1; int N=10;

- The OFELI function banner
outputs the official banner of the library.

banner();

- N is the program argument (if this one is present)
if (argc > 1) N = atoi(argv[1]);

- We now declare an instance of class Mesh
with the appropriate constructor.
Mesh ms(L,N);

- We denote by NbN the number of unknowns, the solution
being prescribed at x=0 and then we print out the mesh. Note that member function
setVerbose with a large value for argument enables printing the whole mesh information.
int NbN = N+1; ms.setVerbose(10); cout << ms;

- We declare an instance of class TrMatrix<double>
for the tridiagonal matrix, with size NbN and an instance of class
Vect<double>
for the right-hand side and the solution. Note that the vector
constructor uses here the mesh instance. This simplifies the use of
algebraic expressions for vector assignments.
TrMatrix<double> A(NbN); Vect<double> b(ms);

- In order to test the code, we choose an exact
solution: We set for right-hand side the function f(x) = 16π
^{2}sin(4πx) and the boundary conditions u(0) = u(1) = 0. This yields the solution u(x) = sin(4πx). To implement this without implementing additional functions we resort to the OFELI's parser. For this, the member function set of class Vect<double> enables assigning regular expressions to nodes in function of their coordinates. The variables are x, y and z.b.set("16*pi*pi*sin(4*pi*x)");

- h is the mesh size (length
of an element).
We can compute the right-hand side of the linear system by multiplying
it by the mesh size.
double h = L/double(N); b *= h;

- We now build up the matrix and the right-hand side. Note that we skip, for the
moment, the first and the last lines for boundary condition treatment. Here x
is the i-th node coordinate.
for (int i=2; i<NbN; i++) { A(i,i ) = 2./h; A(i,i+1) = -1./h; A(i,i-1) = -1./h; }

- We modify the first and last equation in order to take account for boundary conditions.
A(1,1) = 1.; A(1,2) = 0.; b(1) = 0; A(NbN,NbN) = 1.; A(NbN-1,NbN) = 0.; b(NbN) = 0;

- The linear system is solved.
A.solve(b);

- We finally output the solution, calculate the
error at each node and output it.
Note that the exact solution vector use also the parser.
cout << "\nSolution :\n"<< b; Vect<double> sol(ms); sol.set("sin(4*pi*x)"); cout << "Error = " << (b-sol).getNormMax() << endl;

- We now end the program.
return 0; }

- If you execute the program without any argument,
you will obtain as output :
M E S H D A T A =================== Space Dimension : 1 Number of nodes : 11 Number of elements : 10 Number of sides : 0 Solution : 1 0.00000000e+00 2 1.08674760e+00 3 6.71646957e-01 4 -6.71646957e-01 5 -1.08674760e+00 6 -1.62832710e-15 7 1.08674760e+00 8 6.71646957e-01 9 -6.71646957e-01 10 -1.08674760e+00 11 0.00000000e+00 Error = 1.35691088e-01

Copyright © 1998-2016 Rachid Touzani