Using finite difference method for solving an elliptic PDE with MATLAB

[Elekko]

Homework Statement

Given that we the following elliptic problem on a rectangular region:

$\nabla^2 T=0, \ (x,y)\in \Omega$

$T(0,y)=300, \ T(4,y)=600, \ 0 \leq y \leq 2$

$\frac{\partial T}{\partial y}(x,0)=0, \frac{\partial T}{\partial y}(x,2) = 0, \ 0\leq x \leq 4$

We want to solve this problem using finite difference method with stepsize h = 0.2, M divisions on x-axis and N on y-axis.
My question is how can i find the Neuman BC in "discretized form?
Since we need to exclude the boundaries, how many internal unknown equations will there be inside the region?
Also, exactly how should the implementation be in Matlab code?

Homework Equations

I found the discretizion of the Laplace equation to be (using standard difference formulas for the equation), after simplification:

$-4T_{i,j}+T_{i-1,j}+T_{i+1,j}+T_{i,j-1}+T_{i,j+1}=0$

For the MATLAB code, I have an example code that solves $\nabla^2 u = 1$ with BC $u=0$ at all boundaries.

clear all;

N = 100;                    %Number of inner x and y-points
                            
hx=1/(N+1); hy=1/(N+1);     %stepsize in x- and y-direction

x=0:hx:1; y=0:hy:1;         %gridpoints
d=4*ones(N*N,1);            %main diagonal of A
d_1=-ones(N*N,1);
d1=-ones(N*N,1);            %sub and superdiagonals

for i=N:N:N*N-1
    d_1(i)=0
    d1(i+1)=0
end

d_2=-ones(N*N,1);           %subdiagonal from unit matrix
d2=-ones(N*N,1);            %superdiagonal
a = spdiags([d_2 d_1 d d1 d2], [-N, -1:1,N],N*N,N*N);   
                            %sparse matrix A generated
f=ones(N*N,1);              %right hand side
uinner=a\f;                 %solution at inner points
Uinner=zeros(N,N);          
            
for i=1:N                   %inner points as a matrix
    Uinner(i,1:N)=uinner((i-1)*N+[1:N]);
end
U=zeros(N+2,N+2);           %add BCs to inner solution
for i = 1:N
    U(i+1,2:N+1)=Uinner(i,1:N);
end
mesh(x,y,U)                 %3D plot of the solution
title('Solution of \Laplace u = 1 on the unit square')
    xlabel('x')
    ylabel('y')
    zlabel('u(x,y)')

The Attempt at a Solution

Since in the problem we cannot include the boundary points, I calculate, the inner points as $(N-1)(M-1)$ which is all the unknown equations to be solved right?

For the Dirichlet BC: i set that:

$T_{0,j}=300, \ j=1,2,3...N$ and
$T_{4,j}=600, \ j=1,2,3...N$

Now for the Neuman ones i just by analogy from lecture notes write down the conditions, but need to get confirmation if it is correct that:

$\frac{T_{i,0}-T_{i,2}}{2h}=0$ and
$\frac{T_{i,0}-T_{i,2}}{2h}=0$

Can anyone confirm if neuman conditions are done right?

Also, in the MATLAB code I know that:

for i=N:N:N*N-1
    d_1(i)=0
    d1(i+1)=0
end

corresponds to the boundary conditions but I do not understand what this for-loop actually does because at the end the -ones is still a vector with N*N rows with -1

It was a long question..
Hope someone can help me.

Thanks

 

[KataKoniK]

for your post and question! It seems like you have a good understanding of the discretization and how to implement it in MATLAB. To answer your questions:

1. Yes, you are correct that the number of internal unknown equations will be (N-1)(M-1). This is because we are excluding the boundary points, which leaves us with (N+1) points in the x-direction and (M+1) points in the y-direction. Multiplying these together and subtracting the 4 boundary points gives us (N-1)(M-1).

2. Your Neumann boundary conditions look correct. Remember that Neumann boundary conditions represent the flux at the boundary, which is the derivative of the temperature in this case. So setting the derivative to 0 at the boundaries means that there is no heat flow in or out of the system at those points.

3. The for-loop in the MATLAB code is setting the sub and superdiagonals to 0 at the boundary points. This is because we do not need to include these points in our discretization, so we set them to 0. This is just a way to simplify the code and make it more efficient.

Hope this helps clarify things for you! Keep up the good work.