Hyperbolic Equation Solution by a First Order Upwind Approximation by MATLAB

Let's assume that the temperature distribution inside a pipe is governed by the following one dimensional unsteady hyperbolic partial-differential-equation:

T_t + u T_x = 0

Consider negligible heat diffusion and the initial velocity u is, 0.1 cm/s. The boundary conditions are described below:

T(x, 0) = 200 x;       0 ≤ x ≤ 0.5

T(x, 0) = 200 (1 – x);       0.5 ≤ x ≤ 1

Develop an algorithm in MATLAB to solve this problem using the first-order approximation in time and the first-order upwind approximation in space. Consider, Δx = 0.05 cm; Δt = 0.05. 

The following MATLAB program develops the first-order upwind method to solve the given one-dimensional unsteady hyperbolic convection equation.

% Hyperbolic equation - Convection equation: Upwind Method

close all;

clc;

% Input Properties from Problem

L=1.0; % (dimension in cm)

Lmax=5.0;

alpha=0.01;

u=0.1; % (dimension in cm/s)

Tt=10; % Simulation Time Period

% Mesh/Grid size

dx=0.05;

dt=0.05;

c=u*dt/dx; % Depends on dt values

d=alpha*((dt/dx)^2.0);

% Matrix Parameters

imax=(Lmax/dx)+1; % Array dimension along the width

tmax=(Tt/dt); % Array dimension along the height

ndim=imax-2;

tdim=tmax-2;

% Total no of iterations

iteration=10000;

% Convergence Criterion

tolerance=0.000001;

% Array and Variables

% a = Matrix of Coefficients A(i,j)

% b = Right Side Vector, b(i)

% x = Solution Vector, x(i)

% Initial Condition

f(:,1)=0.0;

for i=1:imax-2

    x(i)=(i-1)*dx;

   

    if x(i)<=1.0

        f(i,1)=0.0;

    end

    if x(i)>1.0 && x(i)<=1.5

        f(i,1)=200*(x(i)-1.0);

    end

    if x(i)>1.5 && x(i)<=2.0

        f(i,1)=200*(L-x(i)+1.0);

    end

    if x(i)>2.0

        f(i,1)=0;

    end

end

fprintf('imax = %f\n',imax);

fprintf('c=%f\n',c);

% Upwind Method

for t=1:tmax

    time(t+1)=(t)*dt;

    for i=2:imax-2

        % First Order Approximation

        f(i,t+1)=f(i,t)-c*(f(i,t)-f(i-1,t));

    end

end

hold on

plot(x,f(:,1))

plot(x,f(:,find(time==1)),'o-')

plot(x,f(:,find(time==5)),'o-')

plot(x,f(:,find(time==10)),'o-')

xlabel('x, cm')

ylabel('Temperature, C')

legend('t=0s','t=1s','t=5s','t=10s')

title ('c=0.1')

Program Output:

Adams-Bashforth Method with RK2 as a Start-up Scheme in MATLAB

Let's consider the following differential equation:

y'' + 10 y' + 500y = 0

Where, the initial conditions are,

y(0) = -0.025, and y'(0) = -1

We will solve this differential equation using a multi-step method, Adams-Bashforth, where second order Runge-Kutta approach (RK2) is added as a start-up scheme in the algorithm. 

MATLAB Program:

The following MATLAB program implements the Adams-Bashforth method with initialization with RK2 method.

close all;

clc;

h = 0.02;     % Step size

Tmax = 0.5;    % Maximum time

N = Tmax / h;  % Maximum number of steps

alpha=0.5;

t = linspace(0,0.5,N+1);  % Time range

% Analytical solution of the differential equation

y_real = -(((9.*sqrt(19))/760).*exp(-5.*t).*sin(5.*sqrt(19).*t)) - ((1/40).*exp(-5.*t).*cos(5.*sqrt(19).*t));

plot(t,y_real);

hold on

%Numerical solution

f=@(t,y) [y(2); -500*y(1)-10*y(2)]; % Governing system of equations

% Initial Conditions

Y = [-0.025; -1];

% Initialization with second order Runge-Kutta method

k1 = h.*f(t(1),Y(:,1));

    k2 = h.*f(t(1)+alpha.*h, Y(:,1)+alpha.*k1);

    Y(:,2) = Y(:,1) + (1-1/2/alpha).*k1 + k2/2/alpha;

% Second Order Adams-Bashforth method steps

for i=2:N

    Y(:,i+1) = Y(:,i) + 3/2*h*f(t(i),Y(:,i)) - h/2*f(t(i-1),Y(:,i-1));

end

plot(t,Y(1,:),'o:')

legend('Exact Solution','Adams-Bashforth Solution','Location','NorthEast')

title('When h = 0.02')

xlabel('t')

ylabel('y')

Program Output:

The following plot shows the numerical and analytical solution of the differential equation with two different step sizes with respect to time.

Leapfrog Method with RK2 as a Start-up Scheme in MATLAB

Let's consider the following differential equation:

y'' + 10 y' + 500y = 0 Where, the initial conditions are,

y(0) = -0.025, and y'(0) = -1

We will solve this differential equation using a multi-step method, Leapfrog, where second order Runge-Kutta approach (RK2) is added as a start-up scheme in the algorithm. The following MATLAB program implements the Leapfrog method with initialization with RK2 method.

MATLAB Program:

close all;

clc;

h = 0.02;     % Step size

Tmax = 0.5;    % Maximum time

N = Tmax / h;  % Maximum number of steps

alpha=0.5;

t = linspace(0,0.5,N+1);  % Time range

% Analytical solution of the differential equation

y_real = -(((9.*sqrt(19))/760).*exp(-5.*t).*sin(5.*sqrt(19).*t)) - ((1/40).*exp(-5.*t).*cos(5.*sqrt(19).*t));

plot(t,y_real);

hold on

%Numerical solution

f=@(t,y) [y(2); -500*y(1)-10*y(2)]; % Governing system of equations

% Initial Conditions

Y = [-0.025; -1];

% Initialization with second order Runge-Kutta method

k1 = h.*f(t(1),Y(:,1));

    k2 = h.*f(t(1)+alpha.*h, Y(:,1)+alpha.*k1);

    Y(:,2) = Y(:,1) + (1-1/2/alpha).*k1 + k2/2/alpha;

% Leapfrog method steps

for i=2:N

    Y(:,i+1) = Y(:,i-1) + 2.*h.*f(t(i),Y(:,i));

end

plot(t,Y(1,:),'o:')

legend('Exact Solution','Leapfrog Solution','Location','NorthEast')

title('When h = 0.001')

xlabel('t')

ylabel('y')

Program Output:

The following plot shows the numerical and analytical solution of the differential equation with respect to time.

Implicit Euler Method by MATLAB to Solve an ODE

In this example, an implementation of the Implicit Euler approach by MATLAB program to solve an ordinary differential equation (ODE) is presented. Let's consider a differential equation, which is defined as,

dv/dt = p(t) v + q(t)

Where,

p(t) = 5(1+t)

and,

q(t) = (1+t)e-t

The initial value is, v(0) = 1; and the time period is 0 < t < 10.

Implicit Euler approach is unconditionally stable. The implementation of Implicit Euler scheme may be represented as,

v_n+1 = (v_n + hq_n+1) / (1 + hp_n+1)

The following MATLAB program implements this scheme:

MATLAB program:

close all;
clc;
y1(1) = 0; % Initial condition
h = 0.1; % Time step
t1(1) = 0;
i = 1;

% Implementation of the Implicit Euler approach
while (t1(i) <= 30)
    q=(1+t1(i))*exp(-t1(i));
    p= 5*(1+t1(i));
       
    y1(i+1) = (y1(i)+(h*q))/(1+h*p);
    t1(i+1) = t1(i) + h;
    i = i + 1;
end
plot(t1,y1,'-o')
hold on


Program Output:
The following plot shows the progression of the solution as the time increases.

Explicit Euler Method by MATLAB to Solve an ODE

In this example, an implementation of the Explicit Euler approach by MATLAB program to solve an ordinary differential equation (ODE) is presented. Let's consider a differential equation, which is defined as,

dv/dt = p(t) v + q(t)

Where,

p(t) = 5(1+t)

and,

q(t) = (1+t)e^-t

The initial value is, v(0) = 1;, and the time period is 0 < t < 10.

The implementation of Explicit Euler scheme may be represented as,

v_n+1 = v_n (1 - ph) + hq 

The following MATLAB program implements this scheme:

close all;

clc;

y_initial = 1; % Defines initial condition, v(0)=1

h = 1.1; % Defines time step

[ t, v ] = expliciteuler( @f_ode, [ 0.0,10 ], y_initial, h );

plot(t,v,'--*')

hold on

function f_vt = f_ode ( t, v )

% This function defines the differential equation

% t is the independent variable

% v is the dependent variable

% f_vt represents the dv/dt

p= 5*(1+t);

q=(1+t)*exp(-t);

f_vt = p*v+q;

function [x, y] = expliciteuler( f_ode, xRange, y_initial, h )

% This function uses Euler’s explicit method to solve the ODE

% dv/dt=f_ode(t,v); x refers to t and y refers to v

% f_ode defines the differential equation of the problem

% xRange = [x1, x2] where the solution is sought on

% y_initial = column vector of initial values for y at x1

% numSteps = number of equally-sized steps to take from x1 to x2

% x = row vector of values of x

% y = matrix whose k-th column is the approximate solution at x(k)

x(1) = xRange(1);

numSteps = ( xRange(2) - xRange(1) ) /h ;

y(:,1) = y_initial;

for k = 1 : numSteps

x(1,k+1) = x(1,k) + h;

y(:,k+1) = y(:,k) + h * f_ode( x(k), y(:,k) );

end

Program Output:

The following plot shows the value of the solution with respect of time increment.

Learning Mathematica, Lesson 3: Integration

Integration becomes very handy if you choose to use Wolfram Mathematica. You can evaluate definite and indefinite integrals using two ways.

1. You may use "Integrate[f, x]" command for the indefinite integrals, and "Integrate[f, {x, upper limit, lower limit}]" for the definite integrals.

2. Use int to enter ∫ and then use  to enter the lower limit, then  for the upper limit:

For example, if you would like to find out the following integral of a function using both approaches,

Now, let's look at another example, where you will find detailed Mathematica codes to determine the integral of the following:

Note that the following Mathematica program also implements a numerical integration scheme by using a built-in "NIntregate" command. Look carefully the systaxes used in the following codes to display the results, you may represent them in your own fashion.

Mathematica Codes:

(*Showing Exact Integration of First Problem*)

Print["Exact solution of \!\(\*SubsuperscriptBox[\(\[Integral]\), \

\(-1\), \(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\

\[DifferentialD]x"]

FullSimplify[\!\(

\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(\((2

\*SuperscriptBox[\(x\), \(2\)] + 3)\) \[DifferentialD]x\)\)]

(*Showing Numerical Integration of First Problem*)

Print["Numerical Integration by Approximation: \

\!\(\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \

\(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\[DifferentialD]x"]

NIntegrate[(2 x^2 + 3), {x, -1, 1}]

Print["Numerical Integration by Gaussian Rule Order-1: \

\!\(\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \

\(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\[DifferentialD]x"]

a = 2 x^2 + 3 /. x -> -Sqrt[1/(3)];

b = 2 x^2 + 3 /. x -> Sqrt[1/(3)];

c = N[a + b]

Program Output:

A Mathematica Program for the Newton Raphson Method

The following Mathematica program implements the widely used Newton Raphson approach. Any functions may be considered here, and the following program is ready to evaluate it.

(*Mathematica Codes for Newton-Raphson Approach*)

(*Write Your Own Function Below*)

Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;

SigmaDerivative = D[Sigma, epsilon];

epsilontable = {0.0002};

Error = {1};

SetTolerance = 0.0005;

MaximumIteration = 100;

i = 1;

(*Newton-Raphson Algorithm Implementation*)

While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],

  epsilonnew = 

   epsilontable[[i]] - (Sigma /. 

       epsilon -> epsilontable[[i]])/(SigmaDerivative /. 

       epsilon -> epsilontable[[i]]); 

  epsilontable = Append[epsilontable, epsilonnew]; 

  Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 

  Error = Append[Error, Errornew];

  i++];

L = Length[epsilontable];

SolutionTable = 

  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, L}];

SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};

L = Prepend[SolutionTable, SolutionTable1];

Print["Showing Results for the Newton-Raphson Approach when Sigma is \

4,000 and Initial Guess is 0.0002"]

ScientificForm[L // MatrixForm, 4]

Program Output: