Using finite difference method for solving an elliptic PDE with MATLAB

  • Thread starter Elekko
  • Start date
  • #1
17
0

Homework Statement



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

[itex]\nabla^2 T=0, \ (x,y)\in \Omega[/itex]

[itex]T(0,y)=300, \ T(4,y)=600, \ 0 \leq y \leq 2[/itex]

[itex]\frac{\partial T}{\partial y}(x,0)=0, \frac{\partial T}{\partial y}(x,2) = 0, \ 0\leq x \leq 4[/itex]

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:

[itex]-4T_{i,j}+T_{i-1,j}+T_{i+1,j}+T_{i,j-1}+T_{i,j+1}=0[/itex]

For the matlab code, I have an example code that solves [itex]\nabla^2 u = 1[/itex] with BC [itex]u=0[/itex] at all boundaries.

Code:
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 [itex](N-1)(M-1)[/itex] which is all the unknown equations to be solved right?

For the Dirichlet BC: i set that:

[itex]T_{0,j}=300, \ j=1,2,3...N[/itex] and
[itex]T_{4,j}=600, \ j=1,2,3...N[/itex]

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:

[itex]\frac{T_{i,0}-T_{i,2}}{2h}=0[/itex] and
[itex]\frac{T_{i,0}-T_{i,2}}{2h}=0[/itex]

Can anyone confirm if neuman conditions are done right?

Also, in the matlab code I know that:
Code:
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
 
Last edited:

Answers and Replies

Related Threads on Using finite difference method for solving an elliptic PDE with MATLAB

Replies
0
Views
13K
Replies
1
Views
2K
Replies
5
Views
278
Replies
1
Views
12K
Replies
3
Views
2K
Replies
7
Views
2K
  • Last Post
Replies
3
Views
677
Replies
1
Views
9K
Replies
3
Views
601
Top