Successive over-relaxation of Laplace equation is converging to slightly wrong value. Floating-point error or code error?

39 views (last 30 days)
I am learning how to use the SOR method and I have written a script to find the value of points in a 7x7 grid for the function: psi = sin(x)sinh(y).
The script uses an iterative process using SOR and converges very quickly (as it should for this basic function).
My initial matrix has boundary conditions of sin(x)sinh(1) and sin(1)sinh(y).
HOWEVER it is converging to subtly wrong values (i.e. sin(1/6)sinh(2/6) = 0.05632 and my script is converging to 0.05635).
Why is this happening??
IS THIS A FLOATING-POINT ERROR (if so is there a way to avoid it) OR A CODING ERROR?
Script
gp = 7; %gp^2 is the number of grid points (with the boundary conditions included)
init_psi = zeros(gp,gp);
alpha = 1.35;
N_iter = 30;
for i = 1:gp %defining boundary conditions in a visual way. So first row is at y=1, last row y=0.
init_psi(1,i) = sin((i-1)/(gp-1))*sinh(1);
init_psi(i,end) = sin(1)*sinh(1-((i-1)/(gp-1)));
end
[ psi , hist_values ] = laplace_function(init_psi,alpha,N_iter)
laplace_function
function [ psi , hist_values ] = solve_laplace ( init_psi , alpha , N_iter )
psi = init_psi;
[numRows,numCols] = size(psi); %in this trial it is a square matrix,
%but this function is generalizable for a rectangular grid space
hist_values = zeros(N_iter,3);
lowerRow = round(0.75*numRows); %point in lower half of domain
lowerCol = round(0.25*numCols);
midRow = round(0.5*numRows); %point roughly in middle
midCol = round(0.5*numCols);
upperRow = round(0.25*numRows); %point in upper half of domain
upperCol = round(0.75*numCols);
count = 0;
while count < N_iter %only runs while the number of iterations is less than the maximum
for j = 2:numCols - 1 %loop for 1 iteration
for i = 2:numRows - 1
if j == numCols - 1 %this is necessary so that the code isn't iterating for the defined boundary values
else
Rplus = psi(i,j+1+1)+psi(i,j-1+1)+psi(i-1,j+1)+psi(i+1,j+1)-4*psi(i,j+1);
psi(i,j+1) = psi(i,j+1) + (alpha*Rplus)/4;
end
R = psi(i,j+1)+psi(i,j-1)+psi(i-1,j)+psi(i+1,j)-4*psi(i,j);
psi(i,j) = psi(i,j) + (alpha*R)/4;
end
end %end of iteration
count = count +1;
hist_values(count,1) = psi(lowerRow,lowerCol); %the values in the matrix psi change after each iteration
hist_values(count,2) = psi(midRow,midCol);
hist_values(count,3) = psi(upperRow,upperCol);
end %end of while loop
end %end of function
The output is
which as you can see is slightly off.
Any help would be much appreciated!! :)

Answers (1)

John D'Errico
John D'Errico on 9 Apr 2021
Edited: John D'Errico on 9 Apr 2021
To understand what is happening, you need to first understand what the equations represent!
Is this effectively a partial differential equation on that domain? That is, you are effectively trying to solve an elliptic partial differential equation, the Laplacian, on a simple domain. You have specified boundary conditions that will force the solution to be known.
But are you truly solving the differnetial equation, OR are you making an approximation to the PDE?
It is the latter case. So this is ONLY an approximate solution, in that you have approximated the partial derivatives with finite differnces at each node. Then you solved the resulting linear system of equations for the unknown values at those nodes. (You could also have used backslash here.)
Effectively, it is as if you are solving a simple ODE using Euler's method. What happens when you usea coarse mesh with Euler, versus a fine one? The fine mesh will generate better approximate solution. You need to understand that when you "solve" an ODE or a PDE on a discrete set of points, you are NOT solving the ODE exactly. That is NOT the true solution, but an approximation. You will get a better approximation if you make the mesh finer.
So in reality, this is neither a floating point error, nor is it a coding error. This is an error of approximation, the one thing you did not consider.
And in fact, that is how you would resolve your question. Just use a finer mesh. Does the difference between the ground truth solution and your approximation get smaller? If it does, then you know the answer. This is a typical test for any such numerical approximation.
  1 Comment
Zaphod
Zaphod on 10 Apr 2021
Thanks for your answer! My numerical answer after 30 iterations does indeed get closer to the analytic answer when I reduce the grid spacing. However for a 7x7 grid I was expecting an error of around (as shown here Hatcher, Theodore R. “An Error Bound for Certain Successive Overrelaxation Schemes.” SIAM Journal on Numerical Analysis, vol. 19, no. 5, 1982, pp. 930–941. JSTOR, www.jstor.org/stable/2156985 ). Also a friend of mine achieved this error of using python. So I'm still a bit confused why my error is around .

Sign in to comment.

Categories

Find more on Partial Differential Equation Toolbox in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!