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

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.



Karush–Kuhn–Tucker Conditions (MATLAB)

Code

%KKT with 2 variables, 2 linear constraints with one equality and one

%inequality constraint at once or 1 equality and no inequality constraint at once or no

%equality and one inequality constraint at once

clc;

syms x1 x2 g1 h1 u1 v1 s1

%input of a function like (x1 - 1.5).^2 + (x2 - 1.5).^2 for up to two variables x1, x2

fx=input('Enter the function of ''x'' to be minimized of the form (x - 4).^2 + (y - 6).^2 ');

clc;

C=input('Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint ');

clc;

if C==01

    %constraints input

    g1=input('Enter inequality constraint in the form x1+x2-10 ');

    clc;

    %Lagrange function calculation

    L=fx+(u1*(g1+s1.^2));

    %gradient conditions

    dx1=(diff (L,x1)); dx2=(diff (L,x2)); du1=(diff (L,u1)); ds1=(diff (L,s1));

    %simultaneous solution of gradient conditions

    [sols1,solu1,solx1,solx2]=solve(dx1==0,dx2==0,du1==0,ds1==0);

    a=zeros(3,4);

    for i=1:3

        a(i,1)=solx1(i,1);

        a(i,2)=solx2(i,1);

        a(i,3)=sols1(i,1);

        a(i,4)=solu1(i,1);

    end

    SET1=a(1,:); SET2=a(2,:); SET3=a(3,:);

    %feasibility and non-negativity of Lagrange multiplier check

    if SET1(1,3)>0 || SET1(1,4)>0

        disp('Set 1 '); disp ('     x1    x2    s    u'); disp (SET1)

    end

    if SET2(1,3)>0 || SET2(1,4)>0

        disp('Set 2 '); disp ('     x1    x2    s    u'); disp (SET2)

    end

    if SET3(1,3)>0 || SET3(1,4)>0

        disp('Set 3 '); disp ('     x1    x2    s    u'); disp (SET3)

    end

end

if C==10

    %constraints input

    h1=input('Enter equality constraint in the form x1+x2-10 ');

    clc;

    %Lagrange function calculation

    L=fx+(v1*h1);

    %gradient conditions

    dx1=(diff (L,x1)); dx2=(diff (L,x2)); dv1=(diff (L,v1));

    %simultaneous solution of gradient conditions

    [solv1,solx1,solx2]=solve(dx1==0,dx2==0,dv1==0);

    a=zeros(1,3);

    a(1,1)=solx1;

    a(1,2)=solx2;

    a(1,3)=solv1;

    disp('Set 1 '); disp ('     x1    x2    v'); disp (a);

end

if C==11

    %constraints input

    g1=input('Enter inequality constraint in the form x1+x2-10 ');

    clc;

    h1=input('Enter equality constraint in the form x1+x2-10 ');

    clc;

    %Lagrange function calculation

    L=fx+(v1*h1)+(u1*(g1+s1.^2));

    %gradient conditions

    dx1=(diff (L,x1)); dx2=(diff (L,x2)); du1=(diff (L,u1)); ds1=(diff (L,s1)); dv1=(diff (L,v1));

    %simultaneous solution of gradient conditions

    [sols1,solu1,solv1,solx1,solx2]=solve(dx1==0,dx2==0,du1==0,ds1==0,dv1==0);

    a=zeros(3,5);

    for i=1:3

        a(i,1)=solx1(i,1);

        a(i,2)=solx2(i,1);

        a(i,3)=sols1(i,1);

        a(i,4)=solu1(i,1);

        a(i,5)=solv1(i,1);

    end

    SET1=a(1,:); SET2=a(2,:); SET3=a(3,:);

    %feasibility and non-negativity of Lagrange multiplier check

    if SET1(1,3)>0 || SET1(1,4)>0

        disp('Set 1 '); disp ('       x1         x2         s         u         v'); disp (SET1)

    end

    if SET2(1,3)>0 || SET2(1,4)>0

        disp('Set 2 '); disp ('       x1         x2         s         u         v'); disp (SET2)

    end

    if SET3(1,3)>0 || SET3(1,4)>0

        disp('Set 3 '); disp ('       x1         x2         s         u         v'); disp (SET3)

    end

end

Example 4.29 Arora pg. 128

Step One

In this step the user enters the objective function:

Command Window

Enter the function of 'x' to be minimized (x1 - 1.5).^2 + (x2 - 1.5).^2

Step Two

In this step the user specifies the type(s) and quantity of constraints the objective function is subjected to.

01 means no equality and one inequality constraint, 10 means one equality and no inequality constraint, 11 means one equality and one inequality constraint.

Command Window

Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint 01

Step Three

In this step the user enter the constraint(s).

Command Window

Enter inequality constraint in the form x1+x2-10 x1 + x2 – 2

Step Four

Output of the program; feasible set with candidate minimum point(s).

Command Window

Set 1

     x1    x2    s    u

     1     1     0     1

Example 4.27 Arora pg. 122

Step One

Command Window

Enter the function of 'x' to be minimized (x1 - 1.5).^2 + (x2 - 1.5).^2

Step Two

Command Window

Enter the number of constraints in the problem, e.g. 10 for one equality and no inequality constraint 10

Step Three

Command Window

Enter equality constraint in the form x1+x2-10 x1 + x2 – 2

Step Four

Command Window

Set 1

     x1    x2    v

     1     1     1

Modifications

Extension to five or even more variables and more constraints can be made via simple modifications to the code. The structure of the new code will remain similar to the code mention above, though it will have many more lines of instructions.

LU Decomposition MATLAB

%Please don't mess around with my code
clc;
R=input('How many variables are there in your system of linear equations? ');



clc;
disp ('Please enter the elements of the Coefficient matrix "A" and the Constant/Known matrix "B" starting from first equation');



A=zeros(R,R);

B=zeros(R,1);
for H=1:R

for i=1:R

fprintf ('Coefficient Matrix Column \b'), disp(i), fprintf ('\b\b Row \b'), disp(H);

I=input(' ');



A(H,i)=A(H,i)+I;
end

fprintf ('Known Matrix Row \b'), disp(H);

J=input(' ');



B(H,1)=B(H,1)+J;
end
clc;

[L,U] = lu(A);

Y = L\B;

X = U\Y;
fprintf ('The Coefficient Matrix "A" is \n'), disp(A);

fprintf ('The Constant/Known Matrix "B" is \n'), disp(B);

fprintf ('The Solution Matrix "X" is \n'), disp(X);

Comparison between Down-Force and Drag Produced by a Legacy Spoiler VS a Spoiler with Tubercles (Humpback Whale Fin's Inspired)

Following data was obtained from Simulations carried out in SolidWorks Flow Simulation Premium.

Without Bumps
Air Speed in Km/h	Down Force in N	Drag in N
120	98.682	33.234
110	82.88	27.957
100	68.266	23.02
90	55.299	18.668
80	43.529	14.697
70	33.284	11.255
60	24.438	8.272
50	16.982	5.769
40	10.83	3.688
30	6.08	2.081
20	2.681	0.929
10	0.648	0.235

With Bumps
Air Speed in Km/h	Down Force in N	Drag in N
120	108.238	30.47
110	90.599	25.549
100	74.818	21.047
90	60.423	17.014
80	47.695	13.443
70	36.441	10.27
60	26.682	7.532
50	18.504	5.228
40	11.82	3.352
30	6.613	1.886
20	2.909	0.841
10	0.685	0.211

Comparison between Down Force and Drag

Air Speed in Km/h	Percentage Less Drag		Percentage More Down Force
120	8.32		8.83
110	8.61		8.51
100	8.57		8.76
90	8.86		8.48
80	8.53		8.73
70	8.75		8.66
60	8.95		8.41
50	9.38		8.23
40	9.11		8.38
30	9.37		8.06
20	9.47		7.84
10	10.21		5.4

It is clear that the spoiler with humpback whale's fin's inspired profile not only produce more down force at a particular velocity but also less drag.

Data for Spoiler without Humpback Whale's Fin's Inspired Bumps:

Wing Span: 100 cm

Chord Length: 17.5 cm

Air Velocity: 0-120 Km/h head on

Vertical Pitch: 22.5 Degree Downwards

Gravity Considered

Fluid: Dry Air at STP

Mesh Settings: Coarse (3/10)

Data for Spoiler with Humpback Whale's Fin's Bumps:

Wing Span: 100 cm

Chord Length Large: 17.5 cm

Chord Length Small: 15.75 cm

Air Velocity: 0-120 Km/h head on

Vertical Pitch: 22.5 Degree Downwards

Gravity Considered

Fluid: Dry Air at STP

Mesh Settings: Coarse (3/10)

Let's now take a look at visual representation of data.

This Plot Shows Air Velocity VS Drag, Down-Force by the Spoiler without Bumps

This Plot Shows Air Velocity VS Drag, Down-Force by the Spoiler with Bumps

As you can see from above two plots; the spoiler with the whale's fin like profile generates more down force and less drag.

This Plot Shows Air Velocity VS Down-Force Generated by the Spoilers

The green line represents the Down-Force generated by the spoiler with whale's fin's inspired design. It is around eight percent more at each velocity.

This Plot Shows Air Velocity VS Drag Generated by the Spoilers

The green line represents the Drag generated by the spoiler with whale's fin inspired design. It is around nine percent less at each velocity.

This Plot Shows Air velocity VS Down-Force to Drag Ratio

It is clear from this plot that Down-Force to Drag ratio is around sixteen percent more for whale's fin's inspired spoiler than the legacy one at each velocity.

This Plot Shows Air Flow Around the Spoiler without Bumps at 120 Km/h from the Right Side.

This Plot Shows Air Flow Around the Spoiler without Bumps at 120 Km/h.

This Plot Shows Air Flow Around the Spoiler with bumps at 120 Km/h.

This plot Shows Air Flow Around the Spoiler with bumps at 120 Km/h.

A simple stress analysis was carried out on both spoilers at 120 Km/h. FOS was greater than 1 for both cases.

Advantages of Spoilers:

The main benefit of installing a spoiler on a car is to help it maintain traction at very high speeds. Particularly at speeds around 90 Km/h. A car with a spoiler installed will be easier to handle at highway speeds. Rear spoilers such as the one's analysed in this study; push the back of the car down so the tires can grip the road better and increase stability. It also increases the braking ability of the car.

To build the prototypes and complete the study further, I need donations. To donate your part send an email to fadoobaba@live.com , tweet @fadoobaba, PM at https://www.facebook.com/ThreeDimensionalDesign orhttps://grabcad.com/fahad.rafi.butt or comment with your contact details and I will contact you!. Thank you for reading!

Do comment and share!