Computing the stiffness matrix and the load vector

Computing the stiffness matrix and the load vector

Here are the main commands:

K=Stiffness1(T,fnk)

Assembles the stiffness matrix for the PDE

-div(k*grad u)=f in Omega,

u=0 on Gamma,

du/dn=0 on Bndy(Omega)-Gamma.

The input fnk can be a (positive) number k or a function

implementing a function k(x,y). If fnk is omitted, it is

assumed to be the constant 1. T describes the triangulation

of Omega. For a description of the data structure T, see

"help Mesh1".
help Load1

F=Load1(T,fnf)

Assembles the load vector for the BVP

-div(a*grad u)=f in Omega,

u=0 on Gamma,

du/dn=0 on Bndy(Omega)-Gamma.

The right hand side function f(x,y) must be implemented in the

function fnf.

a=1;

f=@(x,y)x;
Now I compute the stiffness matrix K and the load vector F:
K=Stiffness1(T,a);

F=Load1(T,f);

Finally, I solve the system Ku=F to get the nodal values:

Given the vector of nodal values (and the mesh), you can graph the computed solution using the ShowPWLinFcn1 command:

ShowPWLinFcn1(T,U,g)

This function draws a surface plot of a piecewise linear

function defined on a triangular mesh. The inputs are T,

the mesh (see "help Mesh1" for details about this data structure)

and the vector U, giving the nodal values of the function

(typically U would be obtained by solving a finite element

equation).

The optional input g is a vector of nodal values at the

constrained nodes. If g is not given, this vector is taken

F=Load1(T,f);

u=K\F;

ShowPWLinFcn1(T,u)

K1=Stiffness1(T1,a);

F1=Load1(T1,f);

u1=K1\F1;

There are other utility routines included, such as routines for estimating the and energy norms. For a complete list of routines provided, type "help Fem1".