By Mark S. Gockenbach (siam, 2010)

Testing the code

To see how well the code is working, I can solve a problem with a known solution, and compare the computed solution with the exact solution. I can easily create a problem with a known solution; I just choose a(x,y) and any u(x,y) satisfying the boundary conditions, and then compute

a=@(x,y)1+x.^2;

u=@(x,y)x.*(1-x).*sin(pi*y);
(Notice that u satisfies homogeneous Dirichlet conditions on the unit square.) Then I can compute f:
syms x y

-diff(a(x,y)*diff(u(x,y),x),x)-diff(a(x,y)*diff(u(x,y),y),y)

ans =

2*sin(pi*y)*(x^2 + 1) + 2*x*(x*sin(pi*y) + sin(pi*y)*(x - 1)) - pi^2*x*sin(pi*y)*(x^2 + 1)*(x - 1)

eval(cmd)

f =

@(x,y)2*sin(pi*y)*(x^2+1)+2*x*(x*sin(pi*y)+sin(pi*y)*(x-1))-pi^2*x*sin(pi*y)*(x^2+1)*(x-1)

K=Stiffness1(T,a);

F=Load1(T,f);

U=K\F;

Here is the approximate solution:

ShowPWLinFcn1(T,U)

(The mesh is so coarse that there is only one free node!)

[v,g]=NodalValues1(T,u)

This function sets v equal to the vector of

values of u(x,y) at the free nodes of the mesh T,

and g equal to the vector of values of u(x,y)

at the constrained nodes of the mesh T.

The function implementing u must be vectorized.

See "help Mesh1" for a description of the

data structure for T.
V=NodalValues1(T,u);

Now I will plot the difference between the computed solution and the piecewise linear interpolant of the exact solution (notice the scale on the graph):

K=Stiffness1(T,a);

F=Load1(T,f);

U=K\F;

V=NodalValues1(T,u);

ShowPWLinFcn1(T,V-U)