1D Linear equation by time marching , Finite difference method.
Objective : To solve 1D linear wave equation by time marching method in finite difference using matlab. Observing how the equation diffuses and Analyzing results .
1D linear Wave equation : `(∂u)/(∂t) = c*(∂u)/(∂x)`
Method : Explicit form , using first order forward differencing for time derivative , and 1st oer backward differencing for the space term in the eqn .
Program :
%% 1-D Linear eqn model by time marching method of (Finite difference Method )
% Code by Aadil Shaikh
close all
clear all
clc
%% Specifying parameters
L = 1; % domain length
n = input('Enter number (n) : '); % Numbers
nx = linspace(0,L,n); % Number of steps in grid
T = 0.4 ;
c = 1; % velocity
dt = 1e-2; % time step
dx = nx(2)-nx(1); % width of space
nt = T/dt; % number of steps in time
u = zeros(1,n); % preallocating space
%% define start and end points for square vave
for i = 1:n
if ((0.1<=nx(i))&&(nx(i)<=0.3))
u(i) = 2;
else
u(i) = 1;
end
end
uold = u;
u_initial = u;
%%
for k = 1:nt
%space loop
for i = 2:n
u(i) = uold(i) - (c*dt/dx)*(uold(i)-uold(i-1));
end
%update old velocities
uold = u;
end
%% plotting the velocities
figure(1)
hold all
plot(nx,u_initial,'r','linewidth',1)
plot(nx,u,'-o','linewidth',.5)
axis([0 1.2 0.8 2])
title({['1-D Linear Convection at C = 1 & n = ',num2str(n) ];['time \it = ',num2str(dt*k)]})
xlabel('Grid Co-ordinate (x) \rightarrow')
ylabel('Velocity (u) \rightarrow')
- The parameters are specified as required in the beginning .
- The 'n' (Number) which denotes the number of steps created in the grid between 0-1 , needs to be plotted at different numbers, have been assigned an input to simply ease the process without creating any arrays.
- Basically the results can be plotted now for any number of steps in grid required and is not limited .
- c , time step and dx (width of space) is defined .
- nt is the number of steps in time is a variable utilized to produce the plots at 0.4 seconds .
- a for loop is created to get the velocity profile as the step function , and create the wave , it automatically changes the velocity (u) to 2m/s as required , between the nx value of 0.1 and 0.3
- Another for loop is created , first for loop is the number of steps in time loop and the second one is space loop .
- Where the time marching soln for explicit form is used . using first order forward differencing for time derivative , and 1st oer backward differencing for the space term in the eqn .
- u is the new velocity to be found .
- Plot is created for the initial square wave at the beginning , and then the time marching solution is observed on it . at the req conditions .
- appropriate title and labels are given .
Plot :
n = 20 :
The red square wave is the initial velocity profile .
The blue wave is the final velocity wave profile at time 0.4 seconds and at Space steps number 20
n = 40
The red square wave is the initial velocity profile .
The blue wave is the final velocity wave profile at time 0.4 seconds and at Space steps number 40
n = 80
The red square wave is the initial velocity profile .
The blue wave is the final velocity wave profile at time 0.4 seconds and at Space steps number 80
n = 160
The red square wave is the initial velocity profile .
The blue wave is the final velocity wave profile at time 0.4 seconds and at Space steps number 160 . we can see that wave is distorted and feels like an error . A result cant be obtained from such a wave .
Explanation of plot behaviour :
- As the number of grid spacing is increased one can observe the diffusion of blue wave slows down .
- The reason for square wave not being maintained is , the First order backward differencing scheme used in the space term creates a false diffusion. False diffusion is a type of error observed when the upwind scheme is used to approximate the convection term in convection diffusion equations . False diffusion with the upwind scheme is reduced by increasing the mesh density.
- Basically increasing the number slows down the diffusion . but at higher values it distorts as its result .
- The wave flows to right side , as its being convected by constant linear wave speed c = 1 .
------------------------------------- THE END ------------------------------------------------