# Calculation of error

I have written a code in which I find the approximation of the solution of this elliptic problem.

I calculated the error using the following part of code:

http://pastebin.com/7b5mmuRW

but I get the following errors for $a=-1, b=1, k_1=k_2=1$ :

For N=10:  er = 9.9920e-016
For N=20: er = 2.4425e-015
For N=30: er = 7.1054e-015
For N=40: er = 7.7716e-015
For N=50: er = 1.5765e-014
For N=60: er = 1.6764e-014
For N=70: er = 1.0436e-014


But the error should decrease while $N$ increases. What have I done wrong at the code above?

EDIT:Using this code: http://pastebin.com/crS4vb1t

I get the following results:

[rate,Error]=order

rate =

0.2221   -2.0880   -0.2637   -2.7620


Error =

1.0e-014 *

0.0281    0.0241    0.1024    0.1230    0.8342


EDIT 2:

For $N=4$, $k_1=k_2=1$ the matrix A, with the code that I have written, is the following:

A =

Columns 1 through 14

-4     1     0     0     1     0     0     0     0     0     0     0     0     0
1    -4     1     0     0     1     0     0     0     0     0     0     0     0
0     1    -4     1     0     0     1     0     0     0     0     0     0     0
0     0     1    -4     0     0     0     1     0     0     0     0     0     0
1     0     0     0    -4     1     0     0     1     0     0     0     0     0
0     1     0     0     1    -4     1     0     0     1     0     0     0     0
0     0     1     0     0     1    -4     1     0     0     1     0     0     0
0     0     0     1     0     0     1    -4     0     0     0     1     0     0
0     0     0     0     1     0     0     0    -4     1     0     0     1     0
0     0     0     0     0     1     0     0     1    -4     1     0     0     1
0     0     0     0     0     0     1     0     0     1    -4     1     0     0
0     0     0     0     0     0     0     1     0     0     1    -4     0     0
0     0     0     0     0     0     0     0     1     0     0     0    -4     1
0     0     0     0     0     0     0     0     0     1     0     0     1    -4
0     0     0     0     0     0     0     0     0     0     1     0     0     1
0     0     0     0     0     0     0     0     0     0     0     1     0     0


Columns 15 through 16

 0     0
0     0
0     0
0     0
0     0
0     0
0     0
0     0
0     0
0     0
1     0
0     1
0     0
1     0
-4     1
1    -4


EDIT: Here is the log-log graph of the errors in relation to $N$:

So the errors can't be right. Or am I wrong?

• what you see here is pretty much the machine's error, so it cannot decrease anymore... – Michael Medvinsky Nov 22 '15 at 14:53
• Also, please don't use pastebin for the crucial code, because this question is essentially meaningless once the paste gets deleted (as happened for your other question). – Christian Clason Nov 22 '15 at 14:54
• @MichaelMedvinsky Why does the error decrease for $N=10$ and $N=20$, then for $N=30$ and $N=40$ it increases, then it decreases and increases alternately? – Mary Star Nov 22 '15 at 15:00
• @MichaelMedvinsky Did you see my edit? – Mary Star Nov 22 '15 at 15:27
• @evinda: MichaelMedvinsky already stated the correct answer -- for all practical purposes,1e-15 is zero. So regardless of your $N$, the error is and remains zero. The particular values don't matter -- they're all zero within roundoff. – Wolfgang Bangerth Nov 22 '15 at 15:29

After a many days of discussions, your problems are

1) you choose too simple problem for which you reach the machine zero and therefore cannot observe normal properties 2) perhaps you do not use the boundary values correctly. With the scheme you work the values should be added to the RHS, but since they are zero there is nothing to add. However, one thing to remember, the boundaries points shouldn't be a part of computation.

See the following code, for the two problems (two different sources), one the simple one you worked with and the second one I suggested to solve in comments.

for Kind=1:2 % define what problem to solve
errpre=0;%just to initialize this with something

for N = 1:7

%create grid
n=2^(N+1);
x=linspace(-1,1,n+2);
y=linspace(-1,1,n+2);
h= abs(x(1)-x(2));
[X,Y]=meshgrid(x(2:end-1),y(2:end-1));

%create descrete operator
I = speye(n,n);
E = sparse(2:n,1:n-1,1,n,n);
D = E+E'-2*I;
A = (kron(D,I)+kron(I,D))./(h^2);

%create RHS and the exact solution
if Kind==1
f=2* ((X.^2 - 1) + (Y.^2 - 1));
ex=(X.^2 - 1) .*(Y.^2 - 1);
elseif Kind==2
ex = sin(pi*X).*sin(pi*Y);
f=-ex*2*pi^2;
end
%solve the problem
ap = (A\f(:));

%compute the error
err  = norm(ex(:)-ap(:),inf)/norm(ap(:),inf);

%print the error and the rate
fprintf('Kind=%d\t n=%d\t err=%-10.8d\t rate=%-4.2f\n',Kind,n,err,log2(errpre/err));

%keep the old value for the rate calculation at next step
errpre=err;

end
end


Run this code to get

Kind=1   n=4     err=2.40933816e-16  rate=-Inf
Kind=1   n=8     err=5.69076036e-16  rate=-1.24
Kind=1   n=16    err=8.94357034e-16  rate=-0.65
Kind=1   n=32    err=3.11434148e-15  rate=-1.80
Kind=1   n=64    err=1.11074876e-15  rate=1.49
Kind=1   n=128   err=3.77521199e-15  rate=-1.77
Kind=1   n=256   err=1.42112850e-14  rate=-1.91

Kind=2   n=4     err=1.24859800e-01  rate=-Inf
Kind=2   n=8     err=3.99615153e-02  rate=1.64
Kind=2   n=16    err=1.13319182e-02  rate=1.82
Kind=2   n=32    err=3.01735099e-03  rate=1.91
Kind=2   n=64    err=7.78424525e-04  rate=1.95
Kind=2   n=128   err=1.97680908e-04  rate=1.98
Kind=2   n=256   err=4.98085147e-05  rate=1.99

• I see... Thank you for your answer!!! @MichaelMedvinsky  When we use $k_1 \neq k_2$ do we use the following exact solution and $f$? pastebin.com/M1KziNDC  I get the following results: pastebin.com/kQWQZHyC  Isn't it wrong? – Mary Star Nov 22 '15 at 19:34
• to use $k_1\ne k_2$ you first read carefully my answer to your previous question. As for the exact solution. Unless you have a good book with such problems on your table, it might be difficult to find an exact solution for a singular problem. I would suggest to measure the Cauchy condition of convergence, i.e. take a norm of the difference of the solution on subsequent grids, e.g. $\|u_h-u_{h/2}\|$ where $h$ denote the grid size. – Michael Medvinsky Nov 22 '15 at 19:50
• Since the interval is $[-1,1]$, so it is symmetric, when $N$ is an even number the line $y=0$ isn't on the grid. So we don't have a problem with th interface. Or am I wrong? – Mary Star Nov 22 '15 at 21:15
• yes, but what happens to derivatives across the interface? – Michael Medvinsky Nov 22 '15 at 21:26
• Can we take points for the derivative from the interface although it is not on the grid? – Mary Star Nov 22 '15 at 21:28