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

Question

Fahmid Mahmud 2022-2-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1658800-how-to-calculate-the-errors-for-finite-element-method-for-1d-poisson-equation

评论： Fahmid Mahmud 2022-3-1

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.

2 个评论
显示无隐藏无

Torsten 2022-2-26

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 2022-2-26

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 2022-2-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1658800-how-to-calculate-the-errors-for-finite-element-method-for-1d-poisson-equation#answer_905200

编辑：Torsten 2022-2-26

在 MATLAB Online 中打开

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))

9 个评论
显示 7更早的评论隐藏 7更早的评论

Fahmid Mahmud 2022-2-26

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 2022-2-27

编辑：Torsten 2022-2-27

在 MATLAB Online 中打开

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 2022-3-1

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.

请先登录，再进行评论。

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

2 个评论
显示无隐藏无

回答（1 个）

9 个评论
显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

2 个评论 显示 无隐藏 无

回答（1 个）

9 个评论 显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论