How can solving the boundary value problem using piece-wise linear finite elements on a uniform grid.

I would like to get help or hints to get the codes of Matlab to solve the following quwstion.

[John D'Errico]

@JM

Don't answer your own question just to make a comment about what you did. Moved to a comment:

I tried with the following codes.

***************************************************************************
Variables are:
***************************************************************************
USER INPUTS :
NE         : Number of elements.
NGP        : Number of GQ points used in numerical integration.
aFunc      : a(x) of the model DE.
bFunc      : b(x) of the model DE.
cFunc      : c(x) of the model DE.
fFunc      : f(x) of the model DE.
xMin, xMax : Minimum and maximum coordinates of the problem domain.
BCtype     : BC type of left (min x) and right (max x) boundaries.
Array of size 2.
1: EBC (Dirichlet), 2: NBC (Neumann)
EBCdata    : Array of 2. Stores PV for given EBCs, if there are any.
    NBCdata    : Array of 2. Stores SV for given NBCs, if there are any.
        OTHER VARIABLES :
        NN         : Number of nodes of the mesh. It is NE+1 for linear elements.
        coord      : Coordinates of mesh nodes, array of length NN.
        GQpoint    : Gauss quadrature points, array of size NGP.
        GQweight   : Gauss quadrature weights, array of size NGP.
        ksi        : Coordinate on the reference element.
        S          : Linear shape functions evaluated at GQ points. Array of size
        NGP.
        dS         : Derivatives of shape functions evaluated at GQ points. Array
        of size NGP.
        Ke         : Element level stiffness matrix of size 2x2.
        Fe         : Element level force vector of size 2x1.
        K          : Global stiffness matrix of size NNxNN. Although it is a sparse
        matrix, it is stored in full form for simplicity.
        F          : Global force vector of size NN.
        u          : Array of unknowns. Length NN.
        %}
        clc;
        clear all;
        close all;
        % ************************************************************************
        % Problem dependent data
        % ************************************************************************

Any advice would be appreciated. This is a first time to use Matlab.

[Torsten]

@JM

I think if you use MATLAB for the first time, this task will be too hard to manage.

If you are allowed to use MATLAB code for your problem, use "bvp4c". But I'm in doubt whether it's a finite-element code.

[Bruno Luong]

@Torsten In 1D finite difference and finite element is the same AFAIK.

@JM have you derive the variational or weak form of your equation?

[Torsten]

@Bruno Luong

If in 1d, finite difference and finite element with linear elements is the same, why should the OP derive the variational form of his equations ?

[Bruno Luong]

I guess because it's his home work and it's good to work inside out FEM theory.

Answer 1

[Bruno Luong]

Just for fun (remind me my youth), it actually takes as little as 10 lines of MATLAB to solve it.

% number of dicretization interval or number of elements - 1;
n = 32;
epsilon = 1/10;
% boundary condition values
u0 = 2;
u1 = 2+exp(-1/epsilon);

% Build FEM matrix and rhs
h = 1/n;
S = spdiags(repmat([-1,2,-1]*(epsilon/h)^2 + ...    % S*u = -epsilon^2 u''
                   [1,4,1]*(1/6), n+1, 1), ...      %        +u
                   -1:1, n+1, n+1);
rhs = 1 + h*(0:n)';                                 % 1 + x
% Lock matrix to take into account Dirichlet bc
S(1,2) = 0; rhs(1) = S(1,1)*u0;
S(n+1,n) = 0; rhs(n+1) = S(n+1,n+1)*u1;
% solve for solution
u = S\rhs;

close all
ufun = @(x) 1+x+exp(-x/epsilon);
x = linspace(0,1,n+1);
h = plot(x,u,'.b',x,ufun(x),'r');
legend(h, 'FEM', 'exact')