Different error calculation between Finite Element Method's numerical solution and exact solution

I am comparing the exact and numerical solution (found using Finite Element Method) for sin x over the interval [0, pi]

For n = 8 elements:

numerical solution = [0 0.3728 0.6889 0.9001 0.9742 0.9001 0.6889 0.3728 0];

For n = 16 elements:

numerical solution = [0 0.1938 0.3802 0.5520 0.7026 0.8261 0.9179 0.9745 0.9936 0.9745 0.9179 0.8261 0.7026 0.5520 0.3802 0.1938 0]

I have declared the linspace(0, pi, n+1) to get these values.

Now I want to calculate three error norms using the following formulas:

My exact solution is p(x) = sin (x);

For the first error, my u is a (9x1) vector and my exact solution is a (1x1) sym. So I transformed the exact solution into 9 nodal points by using xvec = (0, pi, 9) and got satisfactory result for the first error:

eps_nodal = norm( (u)' - p(xvec) ) / norm( p(xvec) ) ;

Now I am stuck with the 2nd error, which is the L2 norm. Here I have to transform my u which is a (9x1) vector into a function (my instructor told me piece-wise linear) and get the difference between the function and the sine function. It has to be a function of x because I have to take the derivative of the difference between these 2 function to calculate the 3rd error, H1 norm.

Please let me know how can I convert my u into a function and take the difference between that function and the exact sine function. [****Feel free to ask me if something about the whole question is not clear enough****]

The error norms should have values like this for increasing values of n:

Answer 1

[Torsten]

xvec = linspace(0,pi,9);
f1 = @(x)(interp1(xvec,usol,x,'linear') - sin(x)).^2;
l2_error = sqrt(integral(f1,0,pi))

for the L2-norm.

xvec = linspace(0,pi,9);
dx = xvec(2)-xvec(1);
f2 = @(x)(interp1(xvec,usol,x,'next') - interp1(xvec,usol,x,'previous'))/dx - cos(x)).^2;
h1_error = sqrt(integral(f2,0,pi))

for the H1-norm.

Since there are jumps in the derivative for the calculation of the H1-norm in the discretization points, I'm not sure "integral" will be able to cope with them.

[Fahmid Mahmud]

This seems to be working. Can you please describe in detail what you did here:

f1 = @(x)(interp1(xvec,usol,x,'linear') - sin(x)).^2;

l2_error = sqrt(integral(f1,0,pi))

dx = xvec(2)-xvec(1);

f2 = @(x)(interp1(xvec,usol,x,'next') - interp1(xvec,usol,x,'previous'))/dx - cos(x)).^2;

h1_error = sqrt(integral(f2,0,pi))

[Torsten]

For the l2-norm, the piecewise-linear function you want to use is defined by

@(x)interp1(xvec,usol,x,'linear')

For the h1-norm, the differentiated piecewise-linear function you want to use is defined by

@(x)(interp1(xvec,usol,x,'next') - interp1(xvec,usol,x,'previous'))/dx

If you don't trust the formulae, plot the two functions in [0:pi].