Sunday, 24 April 2022

1D Laplace's equation using Finite Difference Method

This post is about FDM for Laplace Equation with various boundary conditions.

MATLAB Code (1D, Dirichlet Boundary Conditions)

%% initialize the workspace, clear the command window


clear; clc


%% finite difference 1D laplace dirichlet boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, u(L) = 4 %% Ax=b%%


N = 4 ; %number of grid points

L = 1; %length of domain

dx = L/(N-1); %element size


%% initialize variables %%


l = linspace(0,L,N); %independent

u=zeros(1,N); %dependent


%% boundary conditions %%


u(1)=10;

u(end)=4;


%% b vector %%


b=zeros(N-2,1);

b(1) = b(1) - u(1);

b(end) = b(end) - u(end);


%% A matrix


A = -2*eye(N-2,N-2);

for i=1:N-2

    if i<N-2

        A(i,i+1)=1;

    end

    if i>1

        A(i,i-1)=1;

    end

end


%% solve for unknowns %%


x = A\b;


%% fill the u vector with unknowns %%


u(2:end-1) = x;


%% plot the results %%


hold on; grid on; box on, grid minor

set(gca,'FontSize',40)

set(gca, 'FontName', 'Times New Roman')

ylabel('u','FontSize',44)

xlabel('l','FontSize',44)


plot(l,u,'-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)


MATLAB Code (1D, Mixed Boundary Conditions)

%% initialize the workspace, clear the command window

clear; clc

%% finite difference 1D laplace mixed boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, du/dx(L) = 4 %% Ax=b%%

N = 5; %number of grid points
L = 1; %length of domain
dx = L/(N-1); %element size
a = 4;
%% initialize variables %%

l = linspace(0,L,N); %independent
u=zeros(1,N); %dependent

%% dirichlet boundary condition %%

u(1) = 10;

%% b vector %%

b=zeros(N-1,1);
b(1) = b(1) - u(1);
b(end) = b(end) + dx*a; %neumann boundary condition added to b vector

%% A matrix

A = -2*eye(N-1,N-1);
for i=1:N-1
    if i<N-1
        A(i,i+1)=1;
    end
    if i>1
        A(i,i-1)=1;
    end
end
A(N-1,N-1) = -1; %neumann boundary condition added to A matrix
%% solve for unknowns %%

x = A\b;

% fill the u vector with unknowns %%

u(2:end) = x;

%% plot the results %%

hold on; grid on; box on, grid minor
set(gca,'FontSize',40)
set(gca, 'FontName', 'Times New Roman')
ylabel('u','FontSize',44)
xlabel('l','FontSize',44)

plot(l,u','-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)

No comments:

Post a Comment