How to calculate the errors for Finite Element Method for 1D Poisson equation?

I have solved the 1D Poisson equation -u"(x) = sinx using the Finite Element Method over the period 0 to pi.

For n = 8 elements I got this as the (9x1) vector of nodal values:

[0; 0.3827; 0.7071; 0.9239; 1.0; 0.9239; 0.7071; 0.3827; 0]

For n = 16 elements I got this as the (17x1) vector of nodal values:

[0; 0.1951; 0.3827; 0.5556; 0.7071; 0.8315; 0.9239; 0.9808; 1.0; 0.9808; 0.9239; 0.8315; 0.7071; 0.5556; 0.3827; 0.1951; 0]

The Output is like this:

I have used xvec = linspace(0, pi, n+1) for the nodal points (1x9) vector. Now I want to calculate two error norms using the following formulas:

For the double absolute value i.e. norm, I am considering this

double( sqrt(int((error^2), x, 0, pi)) );

The output of errors should be like this in the loglog scale:

Can you help me out with this how to calculate and get the final plot? Thank you.

[Torsten]

Hint:

The exact solution is u*(x) = sin(x).

To get a graph as shown, you will need more than two variations of the number of elements since you can always make a straight line through only two points.

[Fahmid Mahmud]

Thank you for your comment.

I'm just trying to get the error calculation properly with just 2 elements, then I will increase the elements to get the final plot.

Any suggestions on calculating the two errors?

Answer 1

[Torsten]

You must use higher precision for the output of the solution because now, the solution with 16 elements is classified worse than the solution with 8 elements:

u8 = [0; 0.3827; 0.7071; 0.9239;  1.0; 0.9239; 0.7071; 0.3827; 0];
u16 = [0; 0.1951; 0.3827; 0.5556; 0.7071; 0.8315; 0.9239; 0.9808; 1.0; 0.9808; 0.9239; 0.8315; 0.7071; 0.5556; 0.3827; 0.1951; 0];
x8 = linspace(0,pi,numel(u8)).';
x16 = linspace(0,pi,numel(u16)).';
ustar8 = sin(x8);
ustar16 = sin(x16);
eps_nodal8 = norm(u8-ustar8)/norm(ustar8)
eps_nodal16 = norm(u16-ustar16)/norm(ustar16)
eps_L2_8 = trapz(x8,abs(u8-ustar8))
eps_L2_16 = trapz(x16,abs(u16-ustar16))

[Fahmid Mahmud]

I am just not sure how the 16 elements is giving the worse norm. At the same points, the values for 8 and 16 elements are exactly the same, the 16 elements has values in some extra points in the plot. Therefore the 16 elements should have given lower error norm.

Do you have any other suggestion regarding this?

[Torsten]

At the same points, the values for 8 and 16 elements are exactly the same, the 16 elements has values in some extra points in the plot. Therefore the 16 elements should have given lower error norm.

No. The 16 elements vector should give a higher error norm since the error in the remaining 8 element positions is "added" (under the square root).

[Fahmid Mahmud]

You used the norm() function. Can the following be used to get the norm?

double( sqrt(int((error^2), x, 0, pi)) );

And can you explain what your this comment mean? You must use higher precision for the output of the solution because now, the solution with 16 elements is classified worse than the solution with 8 elements

Because solution with 16 elements should be better than 8 elements.

[Torsten]

I made a mistake in the calculation of the L2-error.

The last two lines must read

eps_L2_8 = sqrt(trapz(x8,(u8-ustar8).^2))
eps_L2_16 = sqrt(trapz(x16,(u16-ustar16).^2))

You used the norm() function. Can the following be used to get the norm?

double( sqrt(int((error^2), x, 0, pi)) );

No, because you have no symbolic expressions, but numerical.

And can you explain what your this comment mean? You must use higher precision for the output of the solution because now, the solution with 16 elements is classified worse than the solution with 8 elements

Define

format long

to get the solution with more than 4 decimal places.

[Fahmid Mahmud]

Following your way, I got relatable answers. But I need to calculate a 3rd type of error: H1. The following are all the 3 type of errors that I am considering:

The following is my code for n = 8, 16, 32

n = 8;

y1.u_FEM = [0; 0.3728222934; 0.6888857722; 0.900072637; 0.9742316019; 0.900072637; 0.6888857722; 0.3728222934; 0] % (9x1) vector

xvec1 = linspace(0, pi, 9); % (1x9) vector

error_1 = (y1.u_FEM)' - p(xvec1); % (9x1)' - (1x9) vector

eps_nodal_8 = norm(error_1) / norm(p(xvec1));

eps_L2_8 = sqrt(trapz(xvec1,(error_1).^2));

eps_H1_8 = sqrt(trapz(xvec1,(diff(error_1)).^2));

n = 16;

y2.u_FEM = [0; 0.1938359583; 0.3802229095; 0.5519981074; 0.7025603277; 0.8261235485; 0.9179393047; 0.9744791682; 0.9935703438; 0.9744791682; 0.9179393047; 0.8261235485; 0.7025603277; 0.5519981074; 0.3802229095; 0.1938359583; 0];

xvec2 = linspace(0, pi, 17);

error_2 = (y2.u_FEM)' - p(xvec2);

eps_nodal_16 = norm(error_2) / norm(p(xvec2));

eps_L2_16 = sqrt(trapz(xvec2,(error_2).^2));

eps_H1_16 = sqrt(trapz(xvec2,(diff(error_2)).^2));

n = 32;

y3.u_FEM = [0; 0.0978596622; 0.1947768823; 0.2898182946; 0.3820685982; 0.4706393724; 0.5546776322; 0.6333740432; 0.7059707158; 0.7717685046; 0.8301337408; 0.8805043352; 0.9223951916; 0.955402878; 0.9792095125; 0.993585824; 0.998393361; 0.993585824; 0.9792095125; 0.955402878; 0.9223951916; 0.8805043352; 0.8301337408; 0.7717685046; 0.7059707158; 0.6333740432; 0.5546776322; 0.4706393724; 0.3820685982; 0.2898182946; 0.1947768823; 0.0978596622; 0];

xvec3 = linspace(0, pi, 33);

error_3 = (y3.u_FEM)' - p(xvec3);

eps_nodal_32 = norm(error_3) / norm(p(xvec3));

eps_L2_32 = sqrt(trapz(xvec3,(error_3).^2));

eps_H1_32 = sqrt(trapz(xvec3,(diff(error_3)).^2));

It gives me the following values for the errors:

eps_nodal_error = [0.0257684 0.00642966 0.00160664];

eps_L2_error = [0.0322959 0.00805838 0.00201362];

eps_H1_error = [0 0 0]

Am I doing anything wrong? The H1 error should not be zero. Please advise. Thank you.

[Torsten]

I don't know how d(u_h)/dx is calculated (d(ustar)/dx = cos(x) is clear).

Central differencing, forward differencing, backward differencing - I think it depends on how you approximated u'' = (u')' in the solution of the differential equation.

[Fahmid Mahmud]

Please let me know if the following make sense: For n=8, u is a 9x1 vector consisting of 9 data points. u is used as it is for eps_nodal_error. But for L2 and H1 error, u has to transformed into a piecewise linear function of those 9 data points. I guess then it can be differentiable w.r.t. x {d(u_h)/dx} for H1 error.

Does it make sense?

[Torsten]

The expressions to evaluate must be calculated from the u_h values alone. If you used finite differences instead of finite elements, du_h/dx could be calculated as du_h/dx (xi) ~ (u_h(i+1)-u_h(i-1))/(x(i+1)-x(i-1)) for 2<=i<=n-1, du_h/dx(x1) = (u_h(2)-u_h(1))/(x(2)-x(1)), du_h(xn) = (u_h(n)-u_h(n-1))/(x(n)-x(n-1)). So an expression to compute eps_H1_8 would be

n8 = numel(u8);
x8 = linspace(0,pi,n8).';
du8 = [(u8(2)-u8(1))/(x8(2)-x8(1));(u8(3:n8)-u8(1:n8-2))./(x8(3:n8)-x8(1:n8-2));(u8(n8)-u8(n8-1))/(x8(n8)-x8(n8-1))];
dustar8 = cos(x8);
eps_H1_8 = sqrt(trapz(x8,(du8-dustar8).^2))

But I don't know how the first-order derivatives were approximated within the finite-element method you used. This method should be implemented to calculate du8.

[Fahmid Mahmud]

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.