Search code examples
matlabnumerical-methods

3D Laplace Relaxation in MATLAB

I was able to write a MATLAB program that plots a 1D Laplace relaxation between two metal plates to find equilibrium potential using Jacobi method.

i will be honest, i am not completely sure that i understand what i did, but here is the relevant part of the code:

N=100;
V = zeros(N,N);
V(1,:)=-1;
V(N,:)=1;

for n=1:400
    for i=2:99
        for j=2:99
            V(i,j)=(V(i-1,j)+V(i+1,j)+V(i,j+1)+V(i,j-1))*0.25;
        end
    end
end

and this is how it looks like:

enter image description here

I am wondering if it is possible to do something similar using the same method, but in 3D. I want to visualize something similar in 3D, a 3D potential... box.

I tried using a similar equation that i found in "Computational Physics. N. J. Giordano & H. Nakanishi, Eq.(5.8)" :

 V(i,j,k) = [V(i+1,j,k) + V(i-1,j,k) + V(i,j+1,k) + V(i,j-1,k) + V(i,j,k+1) + V(i,j,k-1)] * (1/6);

and here is the new code that i have been trying to get to work:

N=10; % Used smaller number to reduce processing time.
V = zeros(N,N,N);
V(:,:,1)=-1; %i am using planes instead of axis as "Insulators"
V(:,:,N)=1;

for n=1:100
    for i=2:99
        for j=2:99
            for k=2:99
                V(i,j,k)=(V(i+1,j,k)+V(i-1,j,k)+V(i,j+1,k)+V(i,j-1,k)+V(i,j,k+1)+V(i,j,k-1))*(1/6);
            end
        end
    end
end

And i am getting a Index exceeds matrix dimensions. in the line where V(i,j,k) is.

So again, I am trying to get something that is the 3D version of the 2D graph linked above. *Also, i would really appreciate if someone can explain a little bit (not the MATLAB part, but the math part) what is exactly that i am doing here with that equation, and what this could be used for?

Greatly appreciate the help.

Edit: i forgot to ask: how would you graph this 3D array?

Solution

  • You have defined N but are not using it during the process. Since you have values such as i-1 and j-1 you need to start from 2. And for you have values such as i+1 and j+1 you need to end by N-1. So a working code will look like this:

    N=10; % Used smaller number to reduce processing time.
V = zeros(N,N,N);
V(:,:,1)=-1; %i am using planes instead of axis as "Insulators"
V(:,:,N)=1;

for n=1:1
    for i=2:N-1
        for j=2:N-1
            for k=2:N-1
                V(i,j,k)=(V(i+1,j,k)+V(i-1,j,k)+V(i,j+1,k)+V(i,j-1,k)+V(i,j,k+1)+V(i,j,k-1))*(1/6);
            end
        end
    end
end

    Also, the first for-loop seems to be doing nothing in this piece of code that you have provided.