Solve PDEs using MATLAB (Code)

This is a post about a code I wrote to solve  1-D Poisson's equations. The method employed is called the finite difference method. For 1-D Laplace equation, the RHS of the equation mentioned below goes to zero. Hence all the terms in the B(i) loop equal zero. Also, B(1) and B(Nx) terms equal to the boundary conditions only. I hope this code helps you in learning something.

The example equation is mentioned below.
d²y/dx² = x²+1 and y(0)=10, y(1)=20

The result of the code is a plot between the dependent and the independent variables, as shown in Fig. 1. The accuracy of the solutions lie around 1.5 % of the analytical solution.

%solution to 1D Poisson's equations with dirichlet boundary conditions, 2nd order accurate, uniform grid
clear; clc;
Nx=101; %Enter the number of grid points
clc;
A=zeros(Nx,Nx); %Stiffness matrix
B=zeros(Nx,1); %Dependent variable vector, also used for boundary conditions and other known terms
x=zeros(Nx,1); %Independent variable vector
x(Nx)=1; %Domain size
dx=x(Nx)/(Nx-1); %Grid spacing
B(1)=(((x(1).^2)+1)*(dx.^2))-10; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the first boundary condition
B(Nx)=(((x(Nx).^2)+1)*(dx.^2))-20; %Enter the function the PDE equals to multiplied by the square of the grid spacing minus the second boundary condition
for i=2:Nx-1 %Off boundary values for the dependent variable vector
    B(i)=((x(i).^2)+1)*(dx.^2); %Enter the function the PDE equals to multiplied by the the square of the grid spacing
end
for i=1:Nx-1 %Loop for independent variable
    x(i+1)=(i)*dx;
end
for i=1:Nx %Diagonal terms for the stiffness matrix
    A(i,i)=-2;
end
for i=1:Nx-1 %Off diagonal terms for the stiffness matrix
    A(i,i+1)=1;
    A(i+1,i)=1;
end
X=A\B;
hold on
set(gca,'FontSize',28)
xlabel('Independent Variable','FontSize',38.5)
ylabel('Dependent Variable','FontSize',38.5)
title('Dependent VS Independent Variable (d^2y/dx^2 = x^2+1)','FontSize',38.5);
plot(x,X)

Fig. 1. Plot from the code

Solve ODEs using MATLAB (Code)

I have decided to teach myself some CFD coding using my own intuition and some free online tutorials. This is a post about code I wrote to solve a first order ordinary differential equation. The method employed is called the finite difference method. I will start by coding equations using this method and will slowly move to finite volume method. I hope this code helps you in learning something.

The equation in question is mentioned below.
dy/dx = x+4 and y(0)=3

A uniform grid was employed, as I have not yet figured out how to make non uniform grid. When I figure it out, I will update it here, which will be soon.

The code, mentioned below, starts with a variable called Nx, which represents the total number of grid points in which the domain is divided in to. A rule of thumb, the more points in the domain, more accurate the solution becomes at the cost of computation resources, mainly CPU time and RAM and storage.

x and y are the vectors in which the known and the unknown values are stored. x has the known values, like the length of a pipe through which a fluid flow or a rod through which heat flows etc. y stores the unknown values like temperature of the rod or velocity of the fluid etc. The values are stored at nodes.

y(1) means first value in the vector y. which also happens to be our boundary condition. y is the value of the function at x=0. x=0 is the first value in in the x vector, written as x(1).

x(Nx-1) represents the domain size. Minimum domain size is usually zero e.g. starting point of a rod or a pipe and maximum can be any real number.

To calculate the grid spacing, represented by dx, the domain size was divided by the number of elements. In this example, we have four nodes, two at corners and two in-between, forming a total of three elements.

First order finite difference method was used to discretize the given equation. And a for loop was used for the solution. Using 4 grid points, the solution has an accuracy of 2.22% against the analytical solution of the equation which can be verified by hand calculations. Thank you for reading.

%solution to dy/dx=x+4 y(0)=3
%uniform grid
clear; clc;
Nx=4; %number of grid points
x=zeros(Nx,1); %x vector, known, goes from some known value to maximum value of the domain
y=zeros(Nx,1); %y vector, one value known due to boundary condition rest unknown
y(1)=3; %boundary condition
x(Nx-1)=1; %domain size
dx=x(Nx-1)/(Nx-1); %grid spacing
for i=1:Nx-1
    x(i+1)=(i)*dx;
    y(i+1)=(dx.*(x(i)+4))+y(i); %discretized differential equation
end

Computational Fluid Dynamics: Finite Difference Approximations VS Analytical Differentiation (MATLAB)

Description of the Program

The program takes user input for the function of the variable ‘x’ e.g. x², sin(x), cos(x)/x² etc., the point on which the derivative is evaluated and initial and final grid sizes.

The program then calculates its first derivative analytically.

It then calculates finite difference approximations of the first derivative of order h, with point A unknown and the point A known, using various step sizes, with maximum and minimum step size enter by the user.

 It then calculates and plots the error between the exact value and the approximations of the first derivative.

Code along with Guidelines

clc; %clears the command window

 

%data input and analytical differentiation

 

syms x

fx=input('Enter the function of ''x'' to be differentiated '); %input a function like sin(x) or x.^4 + x.^3

clc; %clears the command window

fxp=(diff (fx,x)); %calculate the derivative analytically

x=input('Enter the point at which the derivative is to be evaluated '); %input the value at which the derivative to be evaluated

clc; %clears the command window

y=inline(fxp,'x'); %convert a statement in to a function so it can be evaluated at some point

z=y(x); %derivative evaluated at the point

 

 

dxmin=input('Enter the smallest value of increment ');

clc;

dxmax=input('Enter the largest value of increment ');

clc;

%finite difference approximation

for dx=dxmin:0.001:dxmax; %loop for different values of grid spacing, dx=h1

    a=1.5;

    h1=a*dx; %h2 in the problem sheet

    h2=(2*a*dx); %h3 in the problem sheet

   

    fc=cos(5+h2)/(5+h2).^2; %value of the function at point c

    fb=cos(5+h1)/(5+h1).^2; %of the function at point b

    fa=cos(5-dx)/(5-dx).^2; %of the function at point a

    fi=cos(x)/(x).^2; %of the function at point A

    fxpu=(fc+fb-(2*fa))/((2*dx)+h1+h2); %formulation with A unknown

    fxpa=(fc+fb+(2*fa)-(4*fi))/((-2*dx)+h1+h2); %formulation with A known

    hold on

    plot ((z-fxpa)/z,dx,'.k') %error plot for approximation with A known

    plot ((z-fxpu)/z,dx,'.r') %error plot for approximation with A unknown

end

%text editing

title('Plot Between Error VS Grid Spacing') %title of the graph

xlabel('Error from the Exact Value'); %x-axis label

ylabel('Step Size/Grid Spacing') %y-axis label

text(0.6, 0.9,'Red: Error plot for approximation with A unknown') %identification

text(0.6, 0.85,'Black: Error plot for approximation with A known') %identification

Calculation

Sample Problem

Command Window after running the code

At Step One

Enter the function of 'x' to be differentiated cos(x)/x.^2

At Step Two

Enter the point at which the derivative is to be evaluated 5

At Step Three

Enter the smallest value of increment 0.001

At Step Four

Enter the largest value of increment 1

Result

 

Observations and Recommendations

It was observed that the error in the finite difference approximation when the value of the function at point A is unknown is less as compared to the case where value of the function at point A is known.

In physical/real-world problems, the exact derivatives of the functions are not known, therefore the error values between the exact and the finite difference approximation cannot be calculated. As it is also observed that the error between finite difference approximation and the exact value of the derivative decreases as the step size decreases, it is advisable to use a step size of ~1/1000 to make a compromise between CPU/RAM usage, time taken to solve the problem and error.