QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP.

54–59, JUL–DEC, 2021 54

An Impact of the Small and Large Grid Sizes on Differential Equations

Abdul Ghafoor Shaikh1,* , Lubna Naz2 , Fahim Raees3 , A.H. Sheikh4 , Asif Ali5
1
Department of Basic Sciences & Related Studies, QUEST Nawabshah
2
Department of Natural Sciences, BNB Women Univesity, Sukkur
3
Mathematics, Faculty of Information Sciences and Humanities, NED UET Karachi
4
Department of Mathematics & Statistics, CCSIS , Institute of Business Management, Karachi
5
Department of Basic Sciences & Related Studies, MUET Jamshoro.
*
Corresponding author: agshaikh@quest.edu.pk

Abstract

The finite difference technique is oldest numerical method to solve differential equations. Like many differential
equations, Helmholtz differential equation which is used to describe many physical phenomena, has long been
solved using finite difference method. can be described by Helmholtz Differential equations. The solution of the
Helmholtz type differential equations is very important. The information that it belongs together because it tells
one coherent story just knowing a little bit about finite differences through to how to solve differential equations an
especial technique is used, how to implement finite difference method and the tool which is used as generic enough
that will immediately be given a whole new differential equation. The analysis of small to moderate sized presented
with the help of a few examples. The improved finite difference method is presented with examples, the method is
simple, clear, and short the MatLab code is available, the improved finite difference method is suitable and easy to
implement, manually as well as computationally.

Keywords—Finite Difference, Differential equations, large grids.

✦

1 Introduction

I n The finite difference approximation method [1],

[2], changes differential equations, whether an ordi-
nary differential or partial differential equation into, a
be shown to implement the finite difference method.
The most important aspect of this paper is that
we will be able to derive a finite difference approxi-
linear system of equations, and it does not give a sym- mation, and know to know how discretization of the
bolic solution. Many problems related to steady state differential equation gives the sparse types matrices
oscillations (mechanical, acoustical, and thermal), and this makes beginners and researchers understand the
wave scattering [3] is modelled by the Helmholtz equa- importance of how, and why we change the governing
tion [4], [5], [6], [7]. Wave scattering has many appli- differential equation into a system of linear equations,
cations in physics, engineering, and science. The basic or sparse matrices. When the talking about finite
philosophy of the finite difference, what would do is in difference approximations that do imply that we don’t
the governing differential equation, we manipulate the have a function a continuous smooth function stored
governing differential equation directly. In the finite in memory, really only know the function at discrete
difference method, we replace the differential notations points, here point f1 , f2 , and f3 , we don’t store the
with some algebraic expressions. The finite difference function between those points and in fact, we don’t
method is a completely numerical method for solv- know what the function is between those points. Some
ing differential equations and working with discrete interpolation to make some good guesses but the point
functions. Finite differences through to how to solve is we don’t know what the function is at in-between
differential equations, using the special technique will points. Suppose we have this discrete function, but
we want to know the slope at some point along this
ISSN: 2523-0379 (Online), ISSN: 1605-8607 (Print) function in this case let’s say want to know the slope
This is an open access article published by Quaid-e-Awam
University of Engineering Science & Technology, Nawabshah,
at f2 . How we do that if only know the function at
Pakistan under CC BY 4.0 International License. discrete points, in this case, will recognize that the
Published on 27.12.2021. derivative is a slope, so we can estimate the at position
QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP. 54–59, JUL–DEC, 2021 55

f2 , by connecting point f1 , and f3 , with a line and we, derivatives at the interior mesh point (i,j) are:
will say that the slope of that line will approximate ∂f fi+1,j − fi−1,j
the slope the true slope at f2 , and the reason we do ( )i,j = + O(h2 )
∂x 2h
that is that our points f1 , and f3 are symmetrically ∂f fi,j+1 − fi,j−1
surrounding the point f2 , and the slope connecting ( )i,j = + O(k 2 )
∂y 2k
points f1 and f3 , which is rise over run how much is our
rise it is f3 − f1 and our run that’s the span between Construct a square matrix so that it premultiplies a
x1 , and x3 , so it’s just 2∆x, where ∆x, is the increment vector, we get a vector containing the first order partial
just from x1 , and x2 . derivative. We would like to figure out what the square
matrix looks when premultiply F, get another column
rise f3 − f1
f ′ (x2 ) = = (1) vector that if we multiply would see the derivative of
run 2∆x F that is, in this case, the first-order derivative.
Similarly the derivatives f ′ (x) for at remaining differ-
ent points x1 , x3 , x4 and x5 will be
rise f2 − f0
f ′ (x1 ) = =
run 2∆x
rise f4 − f2
f ′ (x3 ) = =
run 2∆x
′ rise f5 − f3
f (x4 ) = =
run 2∆x
rise f6 − f4
f ′ (x5 ) = =
run 2∆x
This technique is difficult when considering several
points, which is usual practice to get accurate solution
of the differential equations in that case above method Fig. 2: The derivative of a function goemetricaly pre-
will be cumbersome, to be incorporated, we need to sented at any prescribed point
change strategy;
 
   2 f −f
0
f1 2∆x 
f
  f   3 1 
 −f
   2   2∆x 
     f4 −f2 
?  f3  =  2∆x 
     f5 −f3 
   f4  
2∆x 
f5 f6 −f4
2∆x

The approach by writing a large blank matrix

equation, where our derivative matrix is what we are
calling it, that is what goes in the blank matrix such
that we pre-multiply column vector F we get the
derivative, here we notice that right hand side is the
answer, we have to find out the blank matrix.
 
   2 0 f −f
0 1 0 0 0 f1
−1 0  f2∆x
3 −f1

 1 0 0 f2  
  
2∆x 

   f4 −f2 
 0 −1 0 1 0 f3  =  2∆x  (2)
 
 0

0 −1 0 1 f4  
  f5 −f3 
 2∆x 
0 0 0 −1 0 f5 f6 −f4
2∆x
Fig. 1: Centre Difference Technique Presented at a
point f2 Dealing with a finite difference, we have to keep in
mind and keep control, the first is the positioning of
the points from which the finite-difference approxima-
2 First Order Derivatives tion is calculated, in this case, the more densely they
Using Taylor’s theorem, the second order Central Dif- are packed usually the more accurate, we can estimate
ference Approximations (CDA) to first and second our function or one of its derivatives, but of course
QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP. 54–59, JUL–DEC, 2021 56

that’s more calculation so there is a typical trade-off,

we can improve our accuracy but it requires more cal-
culations which probably requires more memory more
simulation time. In finite difference, the next thing Fig. 3: A 1-D line divided into 10 equal sub-intervals
which is also extremely important is the knowledge
of where we are evaluating our finite difference, and Now different cases are considered with two dif-
we could choose to evaluate the finite difference. It ferent types of boundary conditions and different grid
is also important that, finite differences when we are sizes to compare the solutions at different grid sizes.
using them when we are deriving them we need to
have a good knowledge of the distribution of points 3.1 Case-I
from which we are calculating the finite difference and The grid spacing is denoted by ∆x. The size is taken
where that finite difference is calculating the derivative as N = 11, which is illustrated in Figure 1. For this
or interpolating the function. choice, ∆x can be written as
(xb − xa ) (5 − 0)
∆x = = = 0.5,
3 Problem & Approximating Derivative Us- N −1 (11 − 1)
ing Finite Differences substituting ∆x = 0.5 into the finite difference equa-
Another way to solve the ODE boundary value prob- tion Equation 7 which gives
lems is the finite difference method, where we can use 1 3 2
   
finite difference formulae at evenly spaced grid points + y( i − 1) + 2 − yi
0.52 (2(0.5)) 0.52
to approximate the differential equations. This way, 1 3

we can transform a differential equation into a system +( 2 − y( i + 1) = 0 (8)
0.5 (2(0.5)
of algebraic equations to solve. In the finite difference
and
method, the derivatives in the differential equation
7y( i − 1) − 6yi + y( i + 1) = 0. (9)
are approximated using the finite difference formulae.
We can divide the interval of [a,b] into n equal sub- We have to solve our differential Equation 9, this is a
intervals of length h as presented. simple differential equation, so we have to write, once
We consider a differential equation in every point on the grid, we have eleven points, which
are given
d2 y dy
−3 + 2y = 0 0 ≤ x ≤ 5 (3) 7y3 − 6y4 + y5 = 0 
 
2  y0 = 1
dx dx
7y1 − 6y2 + y3 = 0| 7y4 − 6y5 + y6 (10)
with boundary conditions 
7y2 − 6y3 + y4 = 0 7y5 − 6y6 + y7 = 0


y(0) = 1, y ′ (5) = 5. (4)

 7y6 − 6y7 + y8 = 0 7y9 − 6y10 + y11 + 0 
 

Using the standard finite difference method, we can 7y7 − 6y8 + y9 = 0|

7y8 − 6y9 + y10 = 0 y11 = 5
 
write differential Equation (3) as
(11)
y(x + ∆x) − 2y(x) + y(x − ∆)
and Matrix form will be
∆x2
y(x + ∆x) − y(x − ∆x)
−3 + 2y(x) = 0 (5)
∆x
Conveniently, the finite difference equation can be
written in terms of array indices and this will take the
form
yi+1 − 2yi + yi−1 yi+1−yi−1 2∆x
2
−3 2yi = 0. (6)
∆x +
This makes it easier to implement for computa-
tional purpose;
(12)
1 3 2
   
2
+ y( i − 1) + 2 − yi Solving the matrix equation will give the solution
∆x 2∆x 2∆x2
1 3 (y1 , y2 , y3 , y4 , y5 , y6 , y7 , y8 , y9 , y10 , y11 )

+ − y( i + 1) = 0. (7)
∆x2 2∆x = (1.0, 1.6, 2.5, 3.9, 6.3, 10.0, 15.7, 24.1, 35.1, 41.8, 5.0).
QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP. 54–59, JUL–DEC, 2021 57

4 Improved Finite Difference Scheme The Equation 19 is matrix equation, and boundary
The last one is simpler because we had a simple values are not discretized in the matrix, so the matrix
differential equation it does get more complicated and is simply a raw matrix equation in standard form. We
tedious especially, multiple coupled sets of differential have to incorporate boundary values, the first thing we
equations or pde to solve, so it should be in a manner, want to do is the 1st and last rows, these are the areas
that could implement on MatLab code, working it where we want to force all the values zero, which means
is rather specific to solve that single equation, this we are throwing out the finite difference equation that
new method is much more generic, and can rapidly had written for the 1st and the 11th point, we just
solve different differential equations or different sets of throw them out, and make them all zeros, after that
differential equations, for the small to moderate size we will go into the same rows, and we put a 1 (one)
problem, this is an approved way of implementing the in the diagonal position, then we will go to the last
finite difference method column which had all zeros will insert our boundary
values that is the first point we are forcing value 1, and
y(x) = [y], d/(dx)y(x) = [Dx ], d2 /(dx2 )y(x) = [Dx2 ] then the last point forcing a value of 5.
(13)
Equation 13 takes form;
[Dx2 ][y] − 3[Dx ][y] + 2[y] = 0 (14)
The unknown function of x becomes column vector,
and everything operating y becomes a square matrix,
term by term, and write everything in matrix form,
above Equation 14 is the matrix form of an equation,
but we want to write in terms of a standard form that
is AX = B, for that (20)
([Dx2 ] − 3[Dx ] + 2I)[y] = 0. (15)
Standard form is (y1 , y2 , y3 , y4 , y5 , y6 , y7 , y8 , y9 , y10 , y11 )
[A][y] = 0, (16) = (1.0, 1.6, 2.5, 3.9, 6.3, 10.0, 15.7, 24.1, 35.1, 41.8, 5.0)
where
[A] = ([Dx2 ] − 3[Dx ] + 2I). (17) 4.1 Case-II
This is a way to calculate the square matrix A the sum In this case, the grid size is chosen as ∆x = 0.25.
is an algebraic combination of three matrices, with all Solving the Equation(06), so let us solve this using
same size matrices, what we have to do we will change N = 21 uniform grid points, and spacing ∆x then;
Equation 17 into one line MatLab code i.e.
(xb − xa ) (5 − 0)
[A] = ([Dx2 ] − 3[Dx ] + 2I) = DX2 − 3DX + 2I. (18) ∆x = = = 0.25,
(N − 1) (21 − 1)
This equation is so generic, this becomes very sim- substituting ∆x = 0.25 into the finite difference equa-
ple and rapid to change this matrix for any differential tion Equation(8), which gives
equation, we are solving we would just change it as one
line of code, and that’s neat. Here in this method, we
1 3 2
   
start direct matrix Equation 18, + y( i − 1) + 2 − yi
(0.25) 2 2(0.25) (0.25)2
1 3
 
+ − y( i + 1) = 0 (21)
(0.25)2 2(0.25)

22yi−1 − 30yi + 10yi+1 = 0 (22)

We have to solve the same differential equation with

Equation 21, boundary values so we have simple to
write down the finite difference. Equation 22, once in
(19) every point on the grid, we have twenty one points, this
QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP. 54–59, JUL–DEC, 2021 58

is the sort of raw standard matrix in which boundary

value is not present given.

(26)

4.3 Cale-IV
The different boundary values y(0)= 1 ”and” y(1)=1
are considered in addition to the grid size ∆x = 0.05.
(23)

4.2 Case-III Fig. 5: The 1-D line divided into 20 equal sub-intervals
In this case, the boundary values y(0)=1 and y(5)=2.5
Solving the Equation 06, so let us solve this using
are chosen in addition to the grid size ∆x = 0.125.
N=21 uniform grid points, and spacing ∆x then;
This grid-size is illustrated in Figure 4.
∆x = (xb − xa )/(N − 1) = (1 − 0)/(21 − 1) = 0.05,
substituting ∆x = 0.05 into the finite difference equa-
tion Equation(8), which gives

Fig. 4: The 1-D line divided into 20 equal sub-intervals (1/(0.05)2 + 3/(2.(0.05)))y( i − 1) + (2 − 2/(0.05)2 )yi
+ (1/(0.05)2 − 3/(2.(0.05)))y( i + 1) = 0 (27)

For solving the Equation 06 on grid using N=21;

we have 215y( i − 1) − 399yi + 185y( i + 1) = 0 (28)

∆x = (xb − xa )/(N − 1) = (2.5 − 0)/(21 − 1) = 0.125

Substituting ∆x = 0.125 into the finite difference

equation Equation(8), which gives

(1/(0.125)2 + 3/(2.(0.125)))y( i − 1) + (2 − 2/(0.125)2 )yi

+ (1/(0.125)2 − 3/(2.(0.125)))y( i + 1) = 0 (24)

and

38y( i − 1) − 63yi + 26y( i + 1) = 0 (25) (29)

QUEST RESEARCH JOURNAL, VOL. 19, NO. 2, PP. 54–59, JUL–DEC, 2021 59

[2] J. M. G-Jordan, S. Rojas, M. F-Villegas, and J. E. Castillo,

“A new second order finite difference conservative scheme,”
Divulgaciones Mathemticas, vol. 13, no. 1, pp. 107–122,
2005.
[3] S. Operto, J. Virieux, P. Amestoy, J. L’Excellent, L. Giraud,
and H. Ali, “3D finite-difference frequency-domain model-
ing of visco-acoustic wave propagation using a massively
parallel direct solver: A feasibility study,” Geophysics, vol.
72, no. 5, Art. no. 5, 2007.
[4] A. H. Sheikh, “Development of Helmholtz Solver Based on
Shifted Laplace Preconditioner and a Multigrid Deflation
Technique,” Delft University of Technology, The Nether-
lands, PhD Thesis, 2014.
[5] A. H. Sheikh, C. Vuik, and D. Lahaye, “Fast iterative
solution methods for the Helmholtz equation,” DIAM, TU
Delft, pp. 9–11, 2009.
[6] A. Shaikh G, A. H. Sheikh, A. Asif, and S. Zeb, “Critical
Review of Preconditioners for Helmholtz Equation and their
Spectral Analysis,” Ind. Jour. Sc. Tech, vol. 12, no. 20, May
2019.
[7] A. Bayliss, C. I. Goldstein, and E. Turkel, “On accuracy
conditions for the numerical computation of waves,” Journal
Fig. 6: Comparative between small, moderate and grid of Computational Physics, vol. 59, no. 3, pp. 396–404, 1985.
sizes [8] Y. A. Erlangga and R. Nabben, “On a multilevel Krylov
Method for the Helmholtz Equation preconditioned by
Shifted Laplacian,” Electronic Transactions on Numerical
5 Discussion Analysis (ETNA), vol. 31, pp. 403–424, 2008.
The solutions obtained from different grid sizes as dis- [9] A. H. Sheikh, D. Lahaye, L. Garcia Ramos, R. Nabben, and
C. Vuik, “Accelerating the Shifted Laplace Preconditioner
cussed in different cases in previous section is plotted for the Helmholtz Equation by Multilevel Deflation,” J.
in Figure 3. Curve of solution when smaller grid size is Comput. Phys., vol. 322, no. C, pp. 473–490, Oct. 2016.
chosen is closer to exact solution. Naturally, this would [10] W. A. Siyal, A. H. Sheikh, K. B. Amur, A. Shaikh G, and
R. A. Malookani, “On the Efficiency of Multigrid Solver for
have been the case but at cost of more computational
Shifted Laplace Equation in a Heterogeneous Medium,” Int.
time. For these test problems, the matrix equation J. Appl. Math, vol. 59, no. 3, pp. 102–114, 2020.
is solved directly. However, indirect solver can be
implemented for more large domain problems.

6 Conclusion
The improved finite difference method is tested, which
is specially designed for the hard type of differential
equations, and system of differential equation as well
pde, in matrix oriented form, and how it is built
and assembled, this approach converts the differential
equation into MatLab Codes which makes the hard
type of differential equation(s), easier to be evaluated
as presented. In the last session, the comparison of
different grid sizes and different boundary value prob-
lems was tested and found that small grid sizes are
usually the more accurate, but of course, that needs
more calculation, so there is a typical trade-off, we can
improve our accuracy but it requires more calculations
which probably requires more simulation, more time,
and more memory.

References
[1] I. Singer and E. Turkel, “High-order finite difference meth-
ods for the Helmholtz equation,” Computer Methods in
Applied Mechanics and Engineering, vol. 163, no. 1, Art.
no. 1, 1998.