39 views (last 30 days)

Show older comments

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

John D'Errico
on 9 Apr 2021 at 11:54

Edited: John D'Errico
on 9 Apr 2021 at 11:56

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.

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

Start Hunting!