DISCRETIZATION ERROR ESTIMATION USING THE METHOD OF NEARBY
PROBLEMS: ONE-DIMENSIONAL CASES
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee.
ANIL RAJU
Certificate of Approval:
Anwar Ahmed Christopher. J. Roy, Chair
Associate Professor Assistant Professor
Aerospace Engineering Aerospace Engineering
Brian Thurow Stephen L. McFarland
Assistant Professor Dean
Aerospace Engineering Graduate School
DISCRETIZATION ERROR ESTIMATION USING THE METHOD OF NEARBY
PROBLEMS: ONE-DIMENSIONAL CASES
Anil Raju
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
December 16, 2005.
ii
DISCRETIZATION ERROR ESTIMATION USING THE METHOD OF NEARBY
PROBLEMS: ONE-DIMENSIONAL CASES
Anil Raju
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon request of individuals or institutions and at their expense. The author
reserves all publication rights.
_________________________
Signature of Author
_________________________
Date
iii
VITA
Anil Raju, son of Mr. and Mrs. Thomas Raju, was born on September 3, 1979, in
Trivandrum, India. He graduated with a Bachelors degree in Industrial Engineering from
University of Kerala, Trivandrum, India in May 2002. He joined the graduate program in
Aerospace Engineering at Auburn University in the fall of 2003.
iv
THESIS ABSTRACT
DISCRETIZATION ERROR ESTIMATION USING THE METHOD OF NEARBY
PROBLEMS: ONE-DIMENSIONAL CASES
Anil Raju
Master of Science, December 16, 2005
(Bachelor of Technology, Industrial Engineering, India, May 2002)
91 TYPED PAGES
Directed by Dr. Christopher J. Roy
Discretization error is defined as the difference between the solution of the
discretized equation and the exact solution of the original partial differential equation.
There are two main goals in this study. The first goal is to use of the method of nearby
problems to generate exact solutions to realistic problems so that we can asses the
performance of discretization error estimators can be assessed. The second goal is to
develop and use method of nearby problems itself as an error estimator. Different
polynomial curve fitting techniques are examined and fifth-order Hermite splines are
identified as the best approach for the method of nearby problems. Steady-state Burgers
equation and a modified form of Burgers equation are used as test cases. The analytical
curve fits are then the exact solution to a problem nearby the original problem. Results
are presented for Burgers equation corresponding to a viscous shock wave for Reynolds
v
numbers of 8 and 64, as well as for a modified version of Burgers equation with a
variable viscosity at a nominal Reynolds number of 64. Various discretization error
estimators are evaluated for the original Burgers equation, the nearby problem, and the
modified version of Burgers equation which includes a nonlinear viscosity term. It is
also observed that the method of nearby problems itself performs well as a discretization
error estimator even on coarse meshes.
vi
ACKNOWLEDGEMENT
I would like to thank the Sandia National Laboratories for financial support
extended to this project.
I am sure this is a great opportunity for me to acknowledge several key people I
came across during my educational career in Aerospace Engineering. I would like to
thank my advisor Dr. Christopher J. Roy for honoring me with the project ?Discretization
error estimation using the method of nearby problems: one-dimensional cases?, and
allowing me to work with him. I will be always thankful to him for guiding me and
supervising my work. I will also be obligated to him for the resources he spent on me.
His research and work practices will certainly have significant influence on my career.
I want to express sincere gratitude towards my Masters committee members
Dr.Anwar Ahmed and Dr. Brian Thurow, for their comments and suggestions on my
Masters thesis.
I would like to thank Dr. Matthew Hopkins of Sndia National Laboratories for all
his comments during my work on this project.
Finally, everything I own is dedicated to my Mom, sister Simi, brother-in-law
John, nephew Jeff and last but not the least, Sherin. They mean everything to me.
vii
Anil Raju.
Date: December 16, 2005.
aniraju.c@gmail.com,anilraju2524@yahoo.com
viii
TABLE OF CONTENTS
LIST OF FIGURES ?????????????????????????.. xi
1. INTRODUCTION ????????????????????????? 1
1.1 Background on Computational Fluid Dynamics???????????
1.2 Verification and Validation in CFD????????.. ..??????..
1.3 Sources of Numerical Error??????????????.???...?
1.4 Discretization Error Estimators????????????...????...
1.5 Objective???????.??????????????????....
1
3
4
6
7
2. BACKGROUND?????????????????????????... 8
2.1 Method of Manufactured Solutions???????????.????..
2.2 Prior Work in MNP????.??..?????????????.??
8
10
3. BURGERS EQUATION??????????????????????... 14
3.1 Introduction to Burgers Equation?????....???????????.
3.2 Solutions to Burgers Equation??????????????????.
3.3 Conversion to Dimensional Quantities and Scaling Factors.??????...
14
15
17
4. METHOD OF NEARBY PROBLEMS????????????????? 18
4.1 MNP as an Evaluator of Discretization Error Estimators???????....
4.2 Example for MNP as an Evaluator of Discretization Error Estimators?.......
18
20
5. POLYNOMIAL FITTING PROCEDURES??????????????? 22
5.1 Polynomial Fitting in MNP???????????..???????..
5.2 Standard Polynomial using Matlab?. ??????????????...
5.3 Legendre Polynomial????????????????????..?
5.4 Cubic Spline?????????????????????????
5.5 Fifth Order Hermite Spline??????????????.?????
22
22
24
29
33
6. DISCRETIZATION ERROR ESTIMATORS????????????..?.. 37
ix
6.1 Discretization Error Estimation Using Local Order of Accuracy??..??.
6.2 Discretization Error Estimation Using Global Order of Accuracy?..?..?.
6.3 Mixed Order Error Estimator?????????????????..?
6.4 Method of Nearby Problems???????????????????
37
38
38
39
7. RESULTS????????????????????????????.. 40
7.1 Steady State Burgers Equation??????????????.?..??.
7.2 Nearby Problem to Burgers Equation????????????...?..?.
7.3 Modified Form of Burgers Equation..??????????????..?
7.4 Nearness of the Nearby Problems?????????????????
7.5 Evaluation of Discretization Error Estimates????????????..
40
43
47
48
56
8. CONCLUSIONS?????????????????????????? 68
8.1 Conclusions???????..??????????????.?..??.
8.2 Future Work??????????...???????????...?..?.
68
69
REFERENCES??????????????????????????.?. 70
APPENDICES???????????????????????????? 72
A. Fortran Program for Solving Original Burgers Equation???????..... 72
B. Fortran Program to Solve Nearby Problem to Original Burgers Equation?... 76
C. Fortran Program to Solve Modified Form of Burgers Equation?????? 80
D. Fortran Program to Compute the Coefficients of the Spline Polynomial?.?. 86
E. Matlab Program to Calculate the Source Terms???????????? 90
x
LIST OF FIGURES
1.1: Three Approaches to Fluid Dynamics???????????????.?..
2.1: Norm of source term: Prior work??????.??????????.........
2.2: Norm of source term: Prior work.......?????.????????.???.
3.1: Steady State Exact Solution: Burgers Equation?????????????.
3.2: Unsteady Exact Solution-1: Burgers Equation ???.??????????
3.3: Unsteady Exact Solution-1: Burgers Equation????.?????????.
5.1: Fitting Numerical Solution Using Standard Polynomial: Re=8??????...
5.2: Fitting Numerical Solution Using Standard Polynomial: Re=16??????..
5.3: First Five Legendre Polynomials??.????????????????.
5.4: Fitting Numerical Solution Using Legendre Polynomial: Re=8?..???..??
5.5: Fitting Numerical Solution Using Legendre Polynomial: Re=16??????.
5.6: Fitting Numerical Solution Using Legendre Polynomial: Re=512??????
5.7: Source Term Using Legendre Polynomial Fits, Re=8?..?????????.
5.8: Source Term Using Legendre Polynomial Fits, Re=16.??????????
5.9: Schematic of Spline Fitting System??????????????????
5.10: Fitting Numerical Solution Using Cubic Splines, 9 points: Re=8??????
5.11: Fitting Numerical Solution Using Cubic Splines, 17points: Re=8?????..
5.12: Source Term Using 9 point Cubic Spline Polynomial Fits, Re=8??????
5.13: Source Term Using 17 point Cubic Spline Polynomial Fits, Re=8????..?
5.14: Fitting Numerical Solution Using Hermite Splines, 17 points: Re=8????..
5.15: Fitting Numerical Solution Using Hermite Splines, 65 points: Re=64????
5.16: Fitting Numerical Solution Using Cubic Splines, 129 points: Re=512????
5.17: Source Term Using 9 point Hermite Spline Polynomial Fits, Re=8?????
7.1: Numerical and Exact Solution of Burgers Equation, Re=8?????????
7.2: Numerical and Exact Solution of Burgers Equation, Re=64???????..?
7.3: Discretization Error for Burgers Equation, Re=8???????...?????
7.4: Observed Order of Accuracy: Burgers Equation, Re=8??????????.
7.5: Numerical Solution of Nearby Problem to Burgers Equation, Re=8?????.
7.6: Numerical Solution of Nearby Problem to Burgers Equation, Re=64?????
7.7: Observed Order of Accuracy: Nearby Problem to Burgers Equation, Re=8??..
7.8: Observed Order of Accuracy: Nearby Problem to Burgers Equation, Re=64??
7.9: Solution and Viscosity Variation for Modified Burgers Equation, Re=64???.
7.10: Numerical Solution to Nearby Problem of Modified Burgers Equation???..
7.11: Source Term Using 5 point Hermite Spline Polynomial Fits, Re=8?????.
7.12: Source Term Using 9 point Hermite Spline Polynomial Fits, Re=8?????.
7.13: Source Term Using 17 point Hermite Spline Polynomial Fits, Re=8????...
7.14: Source Term Using 17 point Hermite Spline Polynomial Fits, Re=64????.
7.15: Source Term Using 33 point Hermite Spline Polynomial Fits, Re=64????.
2
12
13
15
16
16
23
23
24
26
26
27
28
28
29
30
31
32
32
34
35
35
36
41
41
42
43
44
45
46
46
47
48
49
49
50
50
51
xi
7.16: Source Term Using 65 point Hermite Spline Polynomial Fits, Re=64????.
7.17: Source Term Using 129 point Hermite Spline Polynomial Fits, Re=512???.
7.18: Source Term Using 257 point Hermite Spline Polynomial Fits, Re=512???.
7.19: Source Term Using 1025 point Hermite Spline Polynomial Fits, Re=512??...
7.20: Source Term of Nearby Problem to Modified Burgers Equation Using 33 point
Hermite Spline Polynomial Fits, Re=64???????????????...
7.21: Source Term of Nearby Problem to Modified Burgers Equation Using 65 point
Hermite Spline Polynomial Fits, Re=64????????????...???
7.22: Source Term of Nearby Problem to Modified Burgers Equation Using 129
point Hermite Spline Polynomial Fits, Re=64??????????..???
7.23: Discretization Error Estimators for Burgers Equation: Re=8, 1025 nodes??..
7.24: Discretization Error Estimators for Burgers Equation: Re=8, 257 nodes??....
7.25: Discretization Error Estimators for Burgers Equation: Re=8, 65 nodes..??..
7.26: Discretization Error Estimators for Burgers Equation: Re=8, 33 nodes???..
7:27: Discretization Error Estimators for Burgers Equation: Re=64, 1025 nodes??
7:28: Discretization Error Estimators for Burgers Equation: Re=64, 257 nodes??..
7:29: Discretization Error Estimators for Burgers Equation: Re=64, 65 nodes???
7:30: Discretization Error Estimators for Burgers Equation: Re=64, 33 nodes???
7.31: Discretization Error Estimators for Nearby Problem, Re=8, 1025 nodes???.
7.32: Discretization Error Estimators for Nearby Problem, Re=8, 257 nodes???..
7.33: Discretization Error estimators for nearby problem, Re=8, 129 nodes????
7.34: Discretization Error Estimators for Nearby Problem, Re=8, 65 nodes????.
7.35: Discretization Error Estimators for Nearby Problem, Re=64, 1025 nodes??..
7.36: Discretization Error Estimators for Nearby Problem, Re=64, 257 nodes???.
7.37: Discretization Error Estimators for Nearby Problem, Re=64, 129 nodes???.
7.38: Discretization Error Estimators for Nearby Problem, Re=64, 65 nodes???..
7.39: Discretization Error Estimators for Nearby Problem to modified Burgers
equation, Re=64, 257 nodes????????????????????..
7.40: Discretization Error Estimators for Nearby Problem to modified Burgers
equation, Re=64, 129 nodes?????????????????..???
51
52
53
53
54
55
55
57
57
58
58
59
60
60
61
62
62
63
63
64
65
65
66
67
67
xii
CHAPTER ONE
INTRODUCTION
1.1 Background on computational fluid dynamics
Computational fluid dynamics (CFD) is a term given to a variety of numerical
techniques applied to solve the equations that govern fluid flow. Before computers,
theory and experiments were the only methods for gaining insight into physical
phenomena. Most of these physical phenomena can be modeled using differential
equations. The ideal approach would be to solve these equations via analytical techniques
and obtain exact solutions. But in most cases, these equations are complex and therefore
difficult to solve analytically. For such cases we have to resort to approximate solutions.
Using computer simulations is one method to obtain these approximate solutions.
The seventeenth century saw the growth of experimental fluid dynamics in
Europe. Experimental fluid dynamics is considered as the first approach [1] in the study
and development of fluid dynamics. The eighteenth and nineteenth centuries saw the
development of the second approach which is theoretical fluid dynamics. The latter part
of the twentieth century saw the development of CFD which is the third approach. The
growth of CFD can be attributed to the development of high-speed computers and
accurate numerical algorithms to solve various problems. CFD complements the other
two approaches of pure theory and pure experiment. There is still a need for theory and
1
experiments, and the future of fluid dynamics depends on the balance of all three
approaches as shown in Fig 1.1.
Pure Pure
Experiment Theory
CFD
Fig 1.1 Three approaches to fluid dynamics
CFD has established itself as a research tool. Moreover, CFD is now establishing
its use as a design tool as well. CFD can be used to predict the presence of vortices in the
flow over vehicles [1] and by studying the behavior of these vortices and their
interaction, one can come up with an optimal aerodynamic design for the vehicle. This is
just one of many uses of CFD as a design tool.
CFD is being used today in a wide variety of areas [1], for example, engine and
automobile applications. By assessing the flow over the body of the vehicle, an
aerodynamic shape can be determined. Today, automotive engineers use CFD to study
details of flow in engines. CFD can be applied in industrial manufacturing as well, for
example, to calculate the flow field in a mold filled by liquid metal or polymer. Civil
engineering uses CFD to tackle problems related to lakes, rivers, estuaries, etc. CFD is
2
also used in other applications relating to heating, air conditioning, and general air
circulation in buildings.
1.2. Verification and validation in CFD
Verification is the process of determining whether the implementation of the
model, accurately represents the original problem of interest [2]. Validation is the process
of determining the extent to which the model accurately represents reality. In other
words, verification checks whether the model is solved correctly while validation ensures
that the correct model is solved. In yet another way, verification can be described as
dealing with the mathematical correctness of the solution, while validation deals with the
physical correctness of the model. Both verification and validation are compared with a
reference standard. In the case of verification the standard is exact solution of the partial
differential equation while in the case of validation, the standard is real world
observations.
Verification is a two step process [3]. The first step is code verification where the
computer code is verified and made free of unacknowledged errors, such as errors due to
mistakes in code or inconsistent numerical schemes. The second step, solution
verification, is the step in which the acknowledged errors such as round-off errors or
discretization errors are assessed.
Code verification [3] can be broadly classified into two areas: numerical
algorithm verification and software quality assurance practices. Some types of algorithm
testing are the method of manufactured solutions, benchmark solutions, iterative
convergence tests, conservation tests, symmetry tests, and so forth. Software quality
3
assurance [4] involves the entire software development process: monitoring the process
and trying to improve it, making sure that all standards and procedures are followed, and
finally making sure that all the defects are found out and properly dealt with. Some of the
software quality analysis tools are static analysis, dynamic testing, and formal testing.
Solution verification has three main aspects [3]. The first aspect, verification of
input data, makes sure that the data that is provided is correct. The second aspect,
numerical error estimation, calculates the acknowledged errors in simulation. The third
aspect, verification of output data, makes sure that the correct post-processing steps are
used.
1.3. Sources of numerical error
In addition to the errors that can come in during the development of the solution
algorithm, there are some acknowledged errors that occur in every computational
simulation. These errors are called numerical errors.
The main types of numerical errors are round-off errors, iterative errors, and
discretization errors. Round-off errors [3] occur due to finite arithmetic in digital
computers. An example of round-off error is ( ) 9999999.00.3/0.10.3 =? for single
precision. Here we see that the computation of in single precision
gives , which leads to the final result of 0.9999999 as opposed to 1.0. Round-
off error can be very important in the case of ill-conditioned problems or time accurate
simulations where a large number of time steps can result in error accumulation. The way
by which the round-off error can be reduced is by using more digits in the computation.
0.3/0.1
3333333.0
4
For example, round-off error can be reduced by using double precision instead of single
precision.
Iterative error [5] is the difference between the current iterative solution and the
exact solution of the discretized equations. When the governing equations for fluid
dynamics are discretized, they often result in a set of non-linear equations. The usual
procedure that is followed to solve these equations is to first linearize them and then
solve them using an iterative method. To stop the iterative process, a convergence
criterion has to be introduced. The iteration goes on until the residual, which is calculated
by substituting the current iterative solution into the discrete equations, is less than this
convergence criterion. If the iterative process is run until the residual is as small as
possible (i.e. machine zero), then the iterative error will be minimized. In this research, a
small iterative error was desired and so a convergence criterion of was used.
14
10
?
Discretization error is defined as the solution of the algebraic system of equations
which is obtained by discretizing the conservation equation and the difference between
the exact solutions of the conservation equations. It is important to estimate the
discretization error before the CFD predictions can be compared with the experimental
data. The first step to solve a set of governing equations numerically is to discretize them.
Discretization is the process of converting the original partial differential equations to an
algebraic set of equations. This algebraic set of equations is then solved on a discrete
mesh to obtain numerical solutions. These solutions are approximate and are generally
different from the exact solution of the governing equations. This difference is the
discretization error. If the discretization approach is consistent, then the discretization
error will decrease as the mesh is refined.
5
1.4. Discretization error estimators
There are a number of ways to estimate the discretization error. Richardson
extrapolation involves the computation of numerical solutions on two or more meshes.
Solutions on these different meshes are then used to compute a higher-order estimate of
the exact solution. This estimate of the exact solution can then be used to estimate the
discretization error. There are certain assumptions that are used in Richardson
extrapolation. The solution is assumed to be smooth, uniform meshes are assumed and
the higher order terms are neglected. The discretization error [3] can be written as
(1.1) TermsOrderHigherhghghgffDE
kkkexactkk
+++=?=
3
3
2
21
where, is the discrete solution on mesh k, is the exact solution to the partial
differential equation, is the coefficient of the i
k
f
exact
f
i
g
th
order term, and h is the measure of the
element size. Consider a second order accurate scheme with solutions on two different
meshes and , with
1
h
2
h
12
2hh = . Neglecting the higher order terms, the discretization
error equation can be written as
(1.2)
2
121
hgff
exact
+=
(1.3)
2
122
)2( hgff
exact
+=
Solving these two equations for and we get
2
g
exact
f
3
21
1
ff
ff
exact
?
+= (1.4)
In general, if we consider a p
th
order accurate scheme with solutions on a fine mesh ( )
and a course mesh ( ), can be approximated as
1
h
2
h
exact
f
6
1
21
1
?
?
+=
pexact
r
ff
ff (1.5)
where r is the grid refinement factor given by
12
/ hhr = . Once is estimated, then the
relative discretization error (RDE) in the fine grid can be calculated as
exact
f
exact
exact
f
ff
RDE
?
=
1
1
(1.6)
1.5. Objective
The current study concentrates on verification. In particular, we focus on solution
verification and use the method of nearby problems (MNP) as an error estimator. We
have also used MNP to come up with an exact solution for a problem that is very close to
(i.e. nearby) the original problem, and have evaluated various error estimators on those
nearby problems.
7
CHAPTER TWO
BACKGROUND
2.1. Method of manufactured solutions (MMS)
MMS [7] provides a general procedure for generating an analytical solution for
code verification. It is a general approach to find coding mistakes/bugs or inconsistent
algorithms. The goal of MMS is code verification. It involves manufacturing an exact
solution to a set of equations which are a modified form of the original partial differential
equations. The solution obtained to this set of modified equations is not physically
realistic. This method is used only to verify the mathematics involved in solving the
original equations, and does not verify the solution obtained by solving the original
equations. This procedure is used when the method of exact solutions cannot be used.
The method of exact solutions [3] is one in which the numerical solution is compared to
an exact solution to the partial differential equation. In the method of exact solution, the
discretization error is computed and then the observed order of accuracy is calculated.
This observed order of accuracy is compared with the formal order of accuracy. This
method [2] is usually not followed for complex cases (geometric complexity, physical
complexity, etc.) because of the limited number of exact solutions.
Once the problem of interest is identified then MMS is conducted for code
verification. MMS is a five step process [3].
8
1. The first step is to choose a manufactured solution. There are certain guidelines
[8] that should be followed in choosing a manufactured solution.
? The manufactured solution should be sufficiently smooth so that the
theoretical order of accuracy can be matched by the observed order of
accuracy on relatively coarse meshes.
? The solution should exercise all the terms of the governing equation. For
example, for the unsteady heat equation, the temperature cannot be chosen
as a function which is independent of time.
? The solution should be such that it has a number of nontrivial derivatives.
For example, in the heat equation which is a second-order equation in
space, picking temperature as a linear function of the spatial coordinate
will not provide a sufficient test.
? The chosen solution should consist of simple analytic functions like
polynomials, trigonometric functions, etc.
? It is better to avoid exponential growth of the solution in time to avoid
confusion with numerical instability.
2. The second step in MMS is to derive the modified governing equation. Here the
governing equations are applied to the chosen manufactured solution. This will
result in the generation of analytical source terms. These analytic source terms
are then added to the governing equations, resulting in a modified form of the
original equations.
3. Once the modified equations are obtained, then they are solved numerically on
different meshes.
9
4. The next step is to evaluate the global discretization error, which can be
calculated as a global norm of the difference between the numerical solution
and the exact solution of the modified equation over the whole domain. The
exact solution of the modified equation is nothing but the manufactured solution
chosen in step 1.
5. The last step in MMS is to compute the observed order of accuracy. If the exact
solution is known, then the observed order of accuracy is calculated as
)ln(
)ln(
1
2
r
DE
DE
p = (2.1)
where, and are the discretization error on the coarse and fine meshes
respectively, and
2
DE
1
DE
r is the grid refinement factor. Once the observed order of
accuracy is calculated, it is compared to the formal order of accuracy. If they
match, then the code is considered to be free of bugs or mistakes in the
discretization. If they do not match, then it suggests that there is some problem
in the code.
2.2. Prior work in MNP
Standard benchmark problems are often used for testing codes. True solutions of
these standard benchmark problems are often known, but it is hard to ascertain the
relationship between the behavior of an algorithm on a standard benchmark problem and
the behavior of the algorithm on the true problem of interest.
Lee and Junkins [9] described the use of a problem near an original ordinary
differential equation (ODE) to serve as a benchmark problem. They constructed a
10
benchmark problem near the original ODE which exactly satisfied the original ODE with
a small, known forcing function. Their work can be summarized as follows:
? Compute a numerical solution on a very refined mesh.
? Generate a polynomial fit for the fine numerical solution by the least squares
approximation using Chebyshev polynomials. A global polynomial was used
instead of a local polynomial to avoid discontinuities.
? Generate the analytic solution from the global fit. This analytic solution
becomes the exact solution of the nearby benchmark problem.
? Use symbolic manipulation to plug the analytic solution into the original
problem and generate small source terms.
? Add these small source terms to the original ODE to form the benchmark
problem.
Since the exact solution of the benchmark problem was known, it was easy to compute
the global error and also to determine optimal integration parameters.
Junkins and Lee [10] later extended the previous methodology to construct exact
special-case solutions for hybrid ODE/PDE systems. This hybrid ODE/PDE system was
able to serve as a benchmark problem to test approximate solution methods. In their
work, they have described a method for coming up with a benchmark problem to
determine optimal time integration parameters, while in our work we use method of
nearby problems to come up with an exact solution to a nearby problem and also as an
error estimator.
Roy and Hopkins [11] examined the generation of exact solutions to problems
near the original problem of interest. They studied two examples: fully developed laminar
11
flow in a channel and a lid-driven cavity. Two codes were used in their work. The
SACCARA code was used to establish a refined numerical solution to the original
problem while the Premo code was used for the implementation of source terms and
generalized boundary conditions. In the first example, the fully developed laminar flow in
a channel they employed a third-order polynomial for the polynomial fit and,
Mathematica was used to compute the analytic source terms. Once the source terms were
computed, they formulated the nearby problem and used the Premo code for the solution
of the nearby problem. Computed solutions on various meshes were compared to the
analytic solution. The L
2
norms of the source term as a function of the underlying mesh
solution are shown in Fig 2.1. It is seen that as the mesh becomes finer (i.e. as h goes to
1), the magnitude of the source terms decrease. Since the magnitude of the source terms
is small, the nearby problem that is formulated will be close to the original problem of
interest. The method of nearby problems was successfully demonstrated for fully
developed laminar flow in a channel.
Fig 2.1. norm of the source term for the third-order polynomial fit (reproduced from
Roy and Hopkins [11])
2
L
12
In the second example, a lid-driven cavity, the analytic source terms were
generated using a fourth-order least square polynomial based on different underlying
mesh refinement levels. To minimize the errors that arise due to singularities at the
corners, a truncated domain was used. The norm of the source term did not get smaller
with mesh refinement as is seen in Fig 2.2. The magnitude of the source term was found
to be large near the boundaries. For this example, the global polynomial did not capture
the solution well.
2
L
Fig 2.2. norm of the source term for the fourth-order polynomial fit (reproduced from
Roy and Hopkins [11])
2
L
13
CHAPTER THREE
BURGERS EQUATION
3.1. Introduction to Burgers equation
Burgers equation is a quasi-linear, one-dimensional, parabolic partial differential
equation of the form
2
2
x
u
x
u
u
t
u
?
?
=
?
?
+
?
?
? (3.1)
where u(x,t) is a two-dimensional scalar field. As an example it can be assumed as
representing velocity as a function of position ?x? and time ?t?, and ??? is the viscosity. A
quasi-linear equation is one in which the highest-order derivative occurs linearly. A
second-order partial differential equation of the form
02
2
22
2
2
=+
?
?
+
?
?
+
?
?
+
??
?
+
?
?
F
y
u
E
x
u
D
y
u
C
yx
u
B
x
u
A (3.2)
is said to be parabolic if the matrix satisfies the determinant
?
?
?
?
?
?
?
C
B
B
A
Z 0=Z
Burgers equation can be related to shock wave theory. The solutions of Burgers
equation can describe the formation and decay of shocks in a compressible fluid. Burgers
equation can also be used as a mathematical model for turbulence. The Navier-Stokes
equation and Burgers equation are quite similar since both contain a non-linear term and
a second-order term multiplied by a small parameter.
3.2. Solution to Burgers equation
14
Benton and Platzmann [12]
describe thirty five solutions to Burgers equation. Out
of the thirty five solutions we have chosen three that are smooth and real:
1. The steady state solution in dimensionless form (denoted by primes) is given by
. This is a solution of Burgers equation when )'tanh(2)','(' xtxu ?= 0
'
=
?
?
t
u
. It
can model a steady shock, and it is smooth, non trivial, and in the real plane.
The plot for this solution is given in Fig. 3.1. for a Reynolds number of 8. The
abscissa is the spatial coordinate x, while the ordinate is the value of u.
Sol 1
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
-5 -4 -3 -2 -1 0 1 2 3 4 5
Sol 1
Fig 3.1. Steady state solution of Burgers equation
2. An unsteady solution to Burgers equation is given as
'
)'cosh(
)'sinh(2
)','('
t
ex
x
txu
?
+
?=
and this solution can model the coalescence of two equal, unsteady shocks and
is shown in Fig. 3.2 for different time levels.
15
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
-6 -4 -2 0 2 4 6
Sol 1.3' t=-6
Sol 1.3' t=-4
Sol 1.3' t=-2
Sol 1.3' t=0
Sol 1.3' t=2
Fig 3.2. Unsteady solution of Burgers equation for shock coalescnese
3. Another unsteady solution to Burgers equation is given as
'4/'2/1
2
'1
'/'
)','('
tx
et
tx
txu
+
= . This solution can model the decay of a solitary pair of
unsteady, equal compression and expansion pulses, as shown in Fig. 3.3.
-4
-3
-2
-1
0
1
2
3
4
-6 -4 -2 0 2 4 6
Sol 2.1(+) t=.1
Sol 2.1(+) t=.2
Sol 2.1(+) t=.5
Sol 2.1(+) t=1
Sol 2.1(+) t=2
Fig 3.3. Unsteady solution of Burgers equation for pulse decay
16
3.3. Conversion to dimensional quantities and scaling factors
The solution of Burgers equation can be converted to dimensional quantities via
transformations given by and . A scaling factor ,/
'
lxx = ,/
2'
ltt ?= ?/
'
ulu = ? can also
be used to scale the solution. The scaling factor was used in the forms ,/?xx =
,/
2
?tt = and uu ?= . The value of the scaling factor depends on the Reynolds number
used. For a Reynolds number of 8, the scaling factor used was 2. The value of the scaling
factor used for a Reynolds number of 64 was 16.
17
CHAPTER FOUR
METHOD OF NEARBY PROBLEMS
The method of nearby problems (MNP) is used to generate exact solutions to
realistic problems, which in turn allows assessment of how discretization error estimators
will perform on the original problem of interest. MNP can also be used as an a posteriori
error estimation technique [13]. It is based on constructing a problem close to the original
problem called the nearby problem. This nearby problem has an exact solution, and also
is representative of the original problem if the source term is sufficiently small. The
nearby problem is numerically solved just like the original problem. Since the exact
solution of the nearby problem is known, the error in its numerical solution can be
evaluated exactly. This information can then be used to estimate the discretization error
in the original problem.
4.1 MNP as an evaluator of discretization error estimators
The method of nearby problems (MNP) involves five steps. These steps can be
summarized as given below:
? Establish an accurate numerical solution
? Generate an analytical curve fit for the above accurate numerical solution
? Generate analytic source terms
? Numerically solve the nearby problem (original problem plus analytical
source term)
18
? Evaluate the discretization error in the nearby problem
An explanation to all of these five steps is given below.
Accurate numerical solution
Once the problem of interest is identified, the first step is to discretize the problem
and come up with an accurate numerical solution.
Analytic curve fit
Once an accurate numerical solution is computed, the second step involves
generating an analytic curve fit to this numerical solution. A curve fitting tool is used to
generate this curve fit. It should be kept in mind that the tool used for the curve fitting
should provide a particular level of continuity which is problem dependent. Once the
curve fit is generated, it should be examined to see how good the fit approximates the
numerical solution. This analytic curve fit will serve as the exact solution to the nearby
problem.
Generation of analytic source terms
The nearby problem differs from the original problem by a (hopefully) small
source term. This source term is obtained by operating the original equation on the
analytic curve fit obtained from the previous step. In the limit as the size of the source
terms approaches zero, the nearby problem approaches the original problem. The
nearness of the nearby problem to the original problem can be judged by calculating the
magnitude of the source term. A more rigorous assessment of the nearness of the nearby
problem is presented in [13] for ordinary differential equations. Such an assessment for
Burgers equation, a nonlinear PDE is beyond the scope of this work.
19
Numerical solution to the nearby problem
The next step involves discretization of the nearby problem and then computation
of the numerical solution on a series of meshes. For a consistent numerical scheme and
sufficiently refined meshes, the formal order of accuracy of the scheme should be
observed, even on the perturbed equations. In general, the discretization error should drop
as , where r is the grid refinement factor and p the formal order of accuracy. In order
to examine the global discretization error behavior, we define the discrete error function
as:
p
r/1
1/2
2
,,
1
1
()
N
k k n exact n
n
E
N
???
=
?
=??
??
?
?
? (4.1)
where k refers to the discrete mesh level and N is the number of mesh nodes in space
including both interior and boundary nodes with the exception of any Dirichlet boundary
nodes for which the discretization error is identically zero. Here, ?
exact,n
refers to the exact
solution (i.e., the curve fit) evaluated at node n.
Evaluation of discretization error
Since the exact solution to the nearby problem is now known, the discretization
error in the numerical solution to the nearby problem no longer has to be estimated, but
can be evaluated exactly.
4.2. Example of MNP as an evaluator of discretization error estimators
The steps involved in MNP can be best explained by simple example. Consider
that
0)(
2
2
=
?
?
+
?
?
=
x
u
x
u
uL (4.2)
20
is the differential operator of interest. The first step in MNP involves obtaining a highly
refined numerical solution to this original problem. Any discretization scheme can be
used, and the numerical solution can be obtained. The second step involves fitting an
analytic curve fit to the refined numerical solution that we have from the first step. This
problem demands continuity as the highest order of the differential operator is two
and the source term that we develop should be slope continuous. So the tool that we use
for curve fitting should be able to provide this continuity criterion. Consider that the
resulting analytic curve fit is
3
C
(4.3)
432
)( exdxcxbxaxu ++++=
?
Now the third step of MNP is to operate the original equation on the analytic curve fit
and come up with an analytic source term. By operating the original problem of interest
on the analytic curve fit, we get
(4.4) )3(4)2(3)1(2)(
2
++++++=
?
xexxdxxcbuL
This is the source term, and it is not equal to zero. Now this becomes the
exact solution of a modified equation or the nearby equation,
)(xs )(
_
xu
)()( xsuL = or 0)()( =? xsuL (4.5)
It should be noted that as s(x) approaches zero, the nearby problem approaches the
original problem. The next step is to come up with a numerical solution to the nearby
problem. Since we have an exact solution to the nearby problem, we can evaluate the
discretization error exactly.
21
CHAPTER FIVE
POLYNOMIAL FITTING PROCEDURES
5.1. Polynomial fitting in MNP
Obtaining a good polynomial fit is the most difficult task in MNP. Two conditions
have to be kept in mind while fitting the numerical solution. The first condition is the
continuity criterion which is problem specific. For our example case, Burgers equation,
continuity is needed in the solution to maintain slope continuity of the source term.
The second condition is that the fit should approximate the numerical solution fairly well
to obtain small source term. The magnitude of the source terms depends on the
approximation that is used. If the approximation that is used is not good, then the
magnitude of the source term will increase and the nearness of the nearby problem is
affected.
3
C
5.2. Standard polynomial using MATLAB
MATLAB [14] uses the function polyfit to fit a polynomial to a given set of data.
Polyfit(X, Y, N) returns the coeffecients of a polynomial P(X) of degree N that fits the
data P(X(I))~=Y(I) in a least-squares sense. Y = Polyval(P,X) gives the value of the
polynomial evaluated at X.
These functions were used to fit a twentieth order polynomial to numerical
solutions of steady-state Burgers equation for various Reynolds numbers. As shown in
Fig 5.1 for the low Reynolds number case, the polynomial fits the data well. But for the
22
high Reynolds number case (Re=64) as shown in Fig 5.2, the global polynomial does not
fit the data well. As a result, the source term will not be sufficiently small, thus this
approach is limited to low Reynolds number (i.e., smoothly varying) cases only.
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Polynomial fit
Numerical solution
Fig 5.1. Fitting the numerical solution of steady state Burgers equation with a standard
20
th
order polynomial: Re=8
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Polynomialfit
Numerical solution
Fig 5.2. Fitting the numerical solution of steady state Burgers equation with a standard
20
th
order polynomial: Re=64
23
5.3. Legendre polynomials
Legendre polynomials are an orthogonal set of polynomials which can be used to
represent a given function. The first few Legendre polynomials are given in equations 5.1
to 5.5:
1)(
0
=xp (5.1)
xxp =)(
1
(5.2)
)13(
2
1
)(
2
2
?= xxp (5.3)
)35(
2
1
)(
3
3
xxxp ?= (5.4)
)33035(
8
1
)(
24
4
+?= xxxp
(5.5)
These first five Legendre polynomials are shown graphically in Fig 5.3.
X
p
i
(x
)
-1 -0.5 0 0.5 1
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Zero
p
0
(x)
p
1
(x)
p
2
(x)
p
3
(x)
p
4
(x)
Fig 5.3. First five Legendre polynomials
24
The Legendre polynomials can be determined by using the following generating function:
()( )
()( )
2/2
0
122!
()
2! ! 2!
k
nkn
n
n
k
nkx
px
knk n k
?
=
??
=
??
?
(5.6)
We can approximate a function f(x) with a truncated Legendre expansion f
n
(x) by
0
() ()
n
ni
i
f xcp
=
=
?
x
(5.7)
where the coefficients c
i
can be found by making use of the orthogonality of the Legendre
polynomials
1
1
1
1
() ()
() ()
i
i
ii
f xp xdx
c
pxpxdx
?
?
=
?
?
(5.8)
The main reason for using Legendre polynomials instead of standard polynomials is that
the Legendre polynomial-based procedure is more stable and robust. That is, the
approximations are guaranteed not to get worse as more terms are included. In Fig 5.4,
the numerical solution of steady-state Burgers equation for Reynolds number of 8 is well
approximated by Legendre polynomials. But when we increase the Reynolds number,
this global polynomial fitting technique again fails to provide a good approximation as is
shown in Fig 5.5 and Fig 5.6. In the presence of a sharp gradient region, the global
polynomial approximation does not perform well.
25
x
u
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical Solution
Legendre Fit
LeastSquare
BurgersEquation
Steady-State
Re= 8
Fig 5.4 Numerical solution and 10
th
order Legendre and standard polynomial fits for
steady-state Burgers equation at Re=8
x
u
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical Solution
Legendre Fit
LeastSquare
BurgersEquation
Steady-State
Re= 16
Fig 5.5 Numerical solution and 10
th
order Legendre and standard polynomial fits for
steady-state Burgers equation at Re=16
26
x
u
-4 -2 0 2 4
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Numerical Solution
Legendre Fit
Leastsquaresfit
BurgersEquation
SteadyState
Re=512
Fig 5.6 Numerical solution and 10
th
order Legendre and standard polynomial fits for
steady-state Burgers equation at Re=512
Since the Legendre polynomial does not approximate the numerical solution properly, the
source term will not be small as shown in the Fig 5.7 for Re=8 and Fig 5.8 for Re=16. In
these cases, the magnitude of the source term is unacceptably large, thus we will not be
able to construct a problem which will be nearby the original problem. So global
polynomials are not good candidates for curve fitting in MNP.
27
x
So
u
r
c
e
T
e
r
m
-4 -3 -2 -1 0 1 2 3 4
-8
-6
-4
-2
0
2
4
6
8
Numerical Solution
10th Order Legendre
8th Order Legendre
6th Order Legendre
4th Order Legendre
BurgersEquation
Steady-State
Re= 8
Fig 5.7 Size of source term using different order Legendre polynomial fits: Re=8
x
So
u
r
c
e
T
e
r
m
-4 -3 -2 -1 0 1 2 3 4
-20
-10
0
10
20
Numerical Solution
10th Order Legendre
8th Order Legendre
6th Order Legendre
4th Order Legendre
BurgersEquation
Steady-State
Re= 16
Fig 5.8 Size of source term using different order Legendre polynomial fits: Re=16
28
5.4. Cubic splines
A cubic spline [15] is a spline polynomial constructed of piecewise third-order
polynomials. A cubic polynomial spline is twice continuously differentiable and depends
on four parameters. It can be written as:
for
32
)()()(:)(
iiiiiiii
xxdxxcxxbaxS ?+?+?+= 1,...,0],,[
1
?=?
+
nixxx
ii
(5.9)
and the setup of the system is given diagrammatically in Fig 5.9. This system has n+1
spline points ranging from 0 to n and n spline zones.
1?n
S
0
S
2?n
S
1
S
0=i ni =1=i 2=i 2?= ni 1?= ni
Fig 5.9. Schematic of the spline fitting system
The conditions that are used to construct the polynomials are
1. niyxS
iii
,...,0,)( ==
2. nixSxS
iiii
,...,1),()(
1
==
?
3. 1,...,1),()(
1
''
?== ? nixSxS
i
i
i
i
4. nixSxS
i
i
i
i ,...,1),()( 1
''''
== ?
The first constraint sets the values at each node. The second constraint sets the continuity
of the values at each node. The third constraint takes care of the continuity of the first
derivative at each interior node and the fourth constraint takes care of the continuity of
29
the second derivative at each interior node. Here we also set
nnn
axS =)( and
for convenience. The first derivatives at end points are also specified
which gives two additional conditions. As we have already seen, global polynomials do
not do a good job approximating sharp gradients. The numerical solution was fit using
cubic splines and the results were encouraging as shown in Figs 5.10, and 5.11.
nn
n cxS 2)(
''
=
Fig 5.10. Fitting the numerical solution with cubic spline fits using 257 nodes and 9
spline points, Re=8
30
Fig 5.11. Fitting numerical solution with cubic spline fits using 257 nodes and 17 spline
points, Re=8
Since the cubic splines perform well in fitting the numerical solution, one can expect that
the source term that result will be small as desired. Plots of the source term using cubic
splines are given in Figs 5.12 and 5.13. Close examination of these figures indicates that
the source terms exhibit slope discontinuities at the spline points. Since cubic splines are
only continuous and Burgers equation contains a second derivative, the first criteria is
not satisfied and we cannot use cubic splines
2
C
31
Fig 5.12. Source term distribution using 9 spline points and 257 nodes
Fig 5.13. Source term distribution using 17 spline points and 257 nodes
32
5.5. Fifth order Hermite spline
A fifth-degree Hermite spline [15] is a spline polynomial constructed of piecewise
fifth-order polynomials. This system also has n+1 spline points ranging from 0 to n and n
spline zones. This can be given by
5432
)()()()()(:)(
iiiiiiiiiiii
xxfxxexxdxxcxxbaxS ?+?+?+?+?+=
for 1,...,0],,[
1
?=?
+
nixxx
ii
(5.10)
where the same spline system shown in Fig. 5.9 is used. The conditions used to construct
a fifth-degree Hermite spline are:
1. niyxS
iii
,...,0,)( ==
2. niyxS i
i
i ,...,0,)(
''
==
3. nixSxS
iiii
,...,1),()(
1
==
?
4. 1,...,1),()( 1
''
?== ? nixSxS
i
i
i
i
5. 1,...,1),()( 1
''''
?== ? nixSxS
i
i
i
i
6. 1,...,1),()( 1
''''''
?== ? nixSxS
i
i
i
i
and also set and as two extra constraints. The first constraint
sets the values at each node. The second constraint sets the first derivative at each node.
The third constraint sets the continuity of the solution at each node. The fourth constraint
takes care of the continuity of the first derivative at each interior node. The fifth
constraint takes care of the continuity of the second derivative at each interior node, and
the sixth constraint sets the third derivative continuity at each interior node.
nnn
axS =)(
nn
n bxS =)(
'
As is seen in Figs 5.14, 5.15 and 5.16, fifth-order Hermite splines do a good job in
approximating the numerical solution for the various Reynolds number cases. Unlike
33
global polynomial fitting techniques, fifth-order Hermite splines approximates even the
high Reynolds number cases well, including sharp gradient region. Fifth-order Hermite
splines were found to be the best fitting tool that can be used for this problem. The
chosen example of Burgers equation demanded continuity which is met by fifth-order
Hermite splines.
3
C
Fig 5.14. Fifth order Hermite spline approximation for Re=8, using 17 spline points
34
Fig 5.15. Fifth order Hermite spline approximation for Re=64, using 65 spline points
Fig 5.16. Fifth order Hermite spline approximation for Re=512, using 129 spline points
35
The numerical solution was fit accurately enough to provide small source terms as
shown in Fig 5.17. for Re=8. In addition, the source term is slope continuous (
continuous) over the entire domain. The size of the source terms for the fifth-order
Hermite splines will be discussed in Chapter 7. A fortran program was developed to
compute the coefficients of the fifth-order Hermite splines (see Appendix D).
1
C
Fig 5.17 Source term distribution using fifth-order Hermite splines with 9 spline
points and 257 nodes, Re=8
36
CHAPTER SIX
DISCRETIZATION ERROR ESTIMATORS
Discretization error arises due to the fact that the original partial differential
equations must be discretized to numerically solve them. There are various methods to
compute this discretization error and in this chapter we discuss those methods.
6.1. Discretization error estimator using local order of accuracy
The discretization error on the fine grid is given as:
exact
ffDE ?=
1
(6.1)
where f
exact
is the exact solution to the partial differential equation and is the
numerical solution. We can estimate f
1
f
exact
using Richardson extrapolation [3] which is
given by:
1
21
1
?
?
+=
p
exact
r
ff
ff (6.2)
where r is the grid refinement factor, p is the order of accuracy, and
2
f and
1
f are the
solutions on the coarse and fine meshes, respectively. Here the order of accuracy can be
computed using solutions on three meshes as
)ln(
ln
12
23
r
ff
ff
p
?
?
?
?
?
?
?
?
?
?
= (6.3)
37
where , and are the solutions on fine, medium and coarse meshes, respectively.
To use Richardson extrapolation with the local order of accuracy as an error estimator,
we need numerical solution on three different meshes.
1
f
2
f
3
f
6.2. Discretization error estimator using global order of accuracy
This technique is different from Richardson extrapolation with local order of
accuracy in the sense that we use the formal order of accuracy instead of computing the
local order of accuracy. As a result, this technique uses numerical solutions on just two
meshes.
6.3. Mixed order error estimator
The mixed order error estimator [16] involves approximation of an exact solution
using three different meshes. In this technique, instead of assuming a single dominant
error term, both first and second order error terms are considered.
(6.4) HOThghgff
kkexactk
+++=
2
21
Three discrete solutions are used to solve a linear set of equations (see [2] for details) and
finally arrive at an approximation for the exact solution given by
2
12
2
23
)1)(1(
))(1()(
?+
??+??
=
rr
ffrrff
f
exact
(6.5)
Once the exact solution is approximated, the discretization error relation is used to
compute the error, which is given as:
exact
ffErrorOrderMixed ?=
1
(6.6)
38
6.4. Method of nearby problems
The Method of nearby problems (MNP) itself can be used as an error estimator.
The discretization error on any given mesh can be evaluated exactly for the nearby
problem as the MNP approach involves the generation of an exact solution to the nearby
problem. If the nearby problem is ?close enough? to the original problem of interest, then
the error on a given mesh for the nearby problem is expected to be very close to the error
in the original problem on the same mesh. The expression for using MNP as an error
estimator is given as:
MNPexactMNP
ffMNP
,,11
?= (6.7)
where is the numerical solution of the nearby problem on a mesh and is
the exact solution of the nearby problem. Here we see that for using this technique, we
use numerical solutions on only one mesh unlike other schemes which use multiple
meshes. MNP thus requires two solutions on the same mesh, one for the original problem
and the second one for the nearby problem.
MNP
f
,1 MNPexact
f
,
39
CHAPTER SEVEN
RESULTS
We have applied MNP to three different one-dimensional test problems:
1. Original Burgers equation (steady state): The exact solution to the steady
state Burgers equation was given by Benton & Platzman [12].
2. ?Nearby? problem to Burgers equation: The exact solution to the nearby
problem to Burgers equation is the fifth order Hermite spline fit to the
numerical solution of the original Burgers equation.
3. Modified form of Burgers equation: A modified form of Burgers
equation was generated which includes a nonlinear viscosity which
varies as a function of both u and x thus giving a nominal Reynolds
number of 64.
7.1. Steady-state Burgers equation
An implicit scheme was used to obtain the numerical solution to Burgers equation
(see Appendix A). Two different Reynolds number cases were run for the original
Burgers equation: Re=8 and Re=64. Both numerical solutions and exact solutions for
these two cases are shown in Figs. 7.1, and 7.2. The numerical solution is right on top of
exact solution which suggests that the implicit scheme used to obtain the numerical
solution is performing well.
40
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical solution
Exact solution
Fig 7.1 Comparison of the numerical solution with exact solution, Re=8
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical solution
Exact solution
Fig 7.2 Comparison of the numerical solution with exact solution, Re=64
41
Since the exact solution of the steady Burgers equation is known, we can evaluate
the discretization error exactly. The discretization error for two different mesh spacing
( and ) is shown in Fig. 7.3. As expected, the discretization error
increases when the mesh spacing increases.
25.0=?x 0625.0=?x
Fig 7.3 Discretization error for the steady state Burgers equation, Re=8
The observed order of accuracy was computed for the various mesh solutions. Since the
exact solution is known, the observed order of accuracy can be computed by the relation
)ln(
)ln(
1
2
r
DE
DE
p = (7.2)
where r is the grid refinement factor and DE
2
and DE
1
are the L
2
norms of the
discretization errors for the coarse and fine mesh, respectively. The plot for the observed
order of accuracy is given in Fig. 7.4. The observed order of accuracy is seen to be
42
approaching two, which is the formal order of accuracy. This suggests that the numerical
solution that was computed is good.
h
O
r
d
e
r
o
f
A
ccu
ra
cy
10
0
10
1
10
2
0
0.5
1
1.5
2
2.5
3
Fig 7.4 Observed order of accuracy for different meshes, Re=8
7.2 Nearby problem to Burgers equation
The nearby problem to Burgers equation is constructed by using fifth-order
Hermite spline to fit the original numerical solution, then operating Burgers equation on
the spline fit to generate analytical source terms. The source terms are then added to the
original problem of interest to give the nearby problem (see Appendix B). When the
source term approaches zero, the nearby problem approaches the original problem. The
same implicit scheme is used to solve the nearby problem. Two different Reynolds
number cases were run for the problem nearby Burgers equation: Re=8 and Re=64. In
both cases an underlying numerical solution from a 1025 node mesh is used. A
43
description of process of choosing the number of spline points is given in detail in
chapter 7.4. The numerical solutions of the nearby problem for Re=8 and Re=64 are
shown in Figs. 7.5 and 7.6. For Re=8 case, 17 spline points were used and for the Re= 64
case, 65 spline points were used. The numerical solution to the nearby problem is
compared with the exact solution to the nearby problem. The figures show that the
numerical solution obtained is good.
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical Solution of nearby problem
Exact solution of nearby problem
Fig 7.5 Numerical solution of the nearby problem to Burgers equation, Re=8
44
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numericalsolutio of nearby problem
Exact solution of nearby problem
Fig 7.6 Numerical solution of the nearby problem to Burgers equation, Re=64
The order of accuracy of the nearby problem was also computed by comparing
numerical solutions on different meshes. The plot of the order of accuracy of the nearby
problem and the original Burger equation is shown in Fig. 7.7 for the Re = 8 case and
Fig. 7.8 for the Re = 64 case. In both the cases, as the mesh becomes fine, the order of
accuracy approaches two which is the formal order of accuracy. This suggests that the
solution that we have obtained is good.
45
x
O
r
d
e
r
o
f
A
ccu
ra
cy
10
0
10
1
10
2
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
2.3
2.4
2.5
OriginalBurgers
Nearby Problem
Fig 7.7 Observed order of accuracy for Burgers equation and the nearby problem, Re=8
x
Or
d
e
r
o
f
A
c
c
u
r
ac
y
10
0
10
1
10
2
1.5
1.6
1.7
1.8
1.9
2
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
3
OriginalBurgers
NearbyProblem
Fig 7.8 Observed order of accuracy for Burgers equation and the nearby problem, Re=64
46
7.3 Modified form of Burgers equation
A modified form of Burgers equation was generated which includes a nonlinear
viscosity which varies as a function of both u and x:
4
1
2
00
4
1
?
?
?
?
?
?
?
?
+
?
?
+
?
?
?
?
?
?
?
?
=
LR
L
xx
xx
u
u
?
?
(7.3)
The constants were chosen as ?
0
= 0.25 m
2
/s, u
0
= 2 m/s, x
L
= - 4 m, and x
R
= 4 m, thus
giving a nominal Reynolds number of 64. This modified form of Burgers equation was
solved numerically using a mesh with 1025 spatial points. The numerical solution and the
viscosity distribution for this modified form of Burgers equation are given in Fig 7.9.
Implicit scheme was used to obtain the numerical solution to the modified form of
Burgers equation (see Appendix C).
x(m)
u(
m
/
s
)
,
?
(m
2
/s
)
-4 -2 0 2 4
-2
-1
0
1
2
u
?
Modified BurgersEquation
Re = 64, Nonlinear Viscosity
Fig 7.9 Solution and viscosity variation for the modified form of Burgers equation
The source terms were computed and added to the modified form of Burgers
equation, which then resulted in a nearby problem the modified form of Burgers
equation. The numerical solution is given in Fig 7.10. The plots shows the numerical
47
solution of modified Burgers equation with the numerical solution of the nearby problem
to modified Burgers equation. The two solutions are very close to each other.
x
So
l
u
t
i
o
n
-4 -2 0 2 4
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Numerical Solution to Nearby Problem of Modfied Burgers
NUmerical Solution to Modified Burgers
Fig 7.10 Numerical solution of the nearby problem to modified form of Burgers equation
7.4. Nearness of the nearby problem
The nearness of the nearby problem can be measured using the source terms.
The source terms are computed by operating the governing equations on the analytical
curve fit. The source term should be small for the nearby problem to be close to the
original problem. As the source term approaches zero, the nearby problem approaches the
original problem of interest. Source terms were calculated for the nearby problem to the
Burgers equation for three different cases: Re=8, Re=64, and Re=512 as shown in Figs.
7.11, 7.12, and 7.13, respectively. In these figures we see that as the number of spline
points increases, the magnitude of the source term decreases. We require that the source
48
term be small. Examining the Figs. 7.11, 7.12, and 7.13, the source term is the smallest
for 17 spline point case and so this case will be the best option.
Fig 7.11 Magnitude of the source term for the nearby problem with 5 spline points, Re=8
Fig 7.12 Magnitude of the source term for the nearby problem with 9 spline points, Re=8
49
Fig 7.13 Magnitude of the source term for the nearby problem with 17 spline points,
Re=8
For the Reynolds number 64 case, the distribution of the source term over the
domain is given in Figs. 7.14, 7.15, and 7.16. Here again as the number of spline points
increases the magnitude of the source term decreases. The smallest source term is
obtained for 65 spline points and so we choose 65 spline points for Re= 64.
Fig 7.14 Magnitude of the source term for the nearby problem with 17 spline points,
Re=64
50
Fig 7.15 Magnitude of the source term for the nearby problem with 33 spline points,
Re=64
Fig 7.16 Magnitude of the source term for the nearby problem with 65 spline points,
Re=64
51
For Re=512, since there is a very sharp gradient region, the number of spline
points to capture the solution should be large. The distribution of the source term for
Re=512 is presented in Figs. 7.17, 7.18, and 7.19. Here we see that for 129 and 257
spline points, the source term is still relatively large and is not going to result in a
problem close to the original problem. But using 1025 spline points, the size of the source
term is decreased. This does not necessarily mean that for higher Reynolds number cases
we have to use a large number of spline points. The number of spline points could be
reduced by the use of spline points with variable spacing. That is, if we can recognize the
areas of large variation, and include a large number of spline points to capture that
variation and fewer points elsewhere, then the number of spline points can be
significantly reduced.
Fig 7.17 Magnitude of the source term for the nearby problem with 129 spline
points, Re=512
52
Fig 7.18 Magnitude of the source term for the nearby problem with 257 spline points,
Re=512
Fig 7.19 Magnitude of the source term for the nearby problem with 1025 spline points,
Re=512
53
Source terms were also computed for the nearby problem to the modified form of
Burgers equation for a Reynolds number of 64. The magnitude of the source terms is
given in Figs. 7.20, 7.21, and 7.22. Here it is found that the magnitude of the source term
decreases when we go from 33 spline point to 65, and then increases slightly when we go
from 65 to 129 spline points. So in this case we should use 65 spline points to formulate
the nearby problem.
Fig 7.20 Magnitude of the source term for the nearby problem to modified Burgers
equation with 33 spline points, Re=64
54
Fig 7.21 Magnitude of the source term for the nearby problem to modified Burgers
equation with 65 spline points, Re=64
Fig 7.22 Magnitude of the source term for the nearby problem to modified Burgers
equation with 129 spline points, Re=64
55
7.5. Evaluation of discretization error estimators
Now that the nearness of the nearby problem has been established, we are in a
position to evaluate the various discretization error estimators. We now examine
discretization error estimates on various grid levels using four different methods: 1)
Richardson extrapolation with global p, i.e. assuming the formal order of accuracy
(requiring two grids), 2) Richardson extrapolation with local p, i.e. Richardson
extrapolation employing the locally calculated order of accuracy (requiring three grids),
3) a mixed-order error estimator (requiring three grids), and 4) the Method of Nearby
Problems (MNP) (requiring only one grid). Numerical solutions are computed on a wide
range of grid levels. In cases where multiple mesh levels are required to obtain the error
estimate, the error is reported for the finest grid only. Grid refinement is performed by
doubling the node spacing (i.e., grid doubling) in all cases.
Fig 7.23 gives the discretization error estimates of all the error estimators
described earlier for meshes using 1025, 513 and 257 points for Re=8 case. For these fine
meshes we see that all the error estimators do a good job in estimating the error. Fig. 7.24
gives the discretization error estimates for meshes using 257, 129 and 65 points. All the
estimators again match the true error. Fig. 7.25 gives the discretization error estimates for
the relatively coarse meshes using 65, 33 and 17 points. The mixed-order error estimator
and Richardson extrapolation using the local order of accuracy underestimates the error,
while MNP and Richardson extrapolation using the global order of accuracy match the
true error. Fig. 7.26 gives the discretization error estimates for very coarse meshes using
33, 17 and 9 points. The mixed-order error estimator and Richardson extrapolation using
local order of accuracy do not provide good estimates, while MNP and Richardson
56
extrapolation using the global order of accuracy are in good agreement with the true
error.
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-3E-06
-2E-06
-1E-06
0
1E-06
2E-06
3E-06
True Error
MNP
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.23 Discretization error estimators for Burgers equation, Re=8, finest mesh =1025
points
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-5E-05
-4E-05
-3E-05
-2E-05
-1E-05
0
1E-05
2E-05
3E-05
4E-05
5E-05
True Error
MNP
MixedOrder
RDE withlocalp
RDE withglobalp
Fig 7.24 Discretization error estimators for Burgers equation, Re=8, finest mesh = 257
points
57
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-0.0008
-0.0006
-0.0004
-0.0002
0
0.0002
0.0004
0.0006
0.0008
True Error
MNP
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.25 Discretization error estimators for Burgers equation, Re=8, finest mesh = 65
points
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-0.003
-0.002
-0.001
0
0.001
0.002
0.003
True Error
MNP
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.26 Discretization error estimators for Burgers equation, Re=8, finest mesh = 33
points
58
Next we examine the higher Reynolds number case with Re=64. Here also we
expect that all the error estimators will perform well on the finer meshes, while on
coarser meshes, some of the error estimators will break down. The plots of the
discretization error estimators for the Reynolds number of 64 are given in Figs. 7.27,
7.28, 7.29, and 7.30. In Fig 7.27, where the finest mesh used is 1025 points, all the error
estimators perform well. In Fig 7.28, where the finest mesh is 257 points, the mixed-order
error estimator and Richardson extrapolation with the global order of accuracy begin to
deviate from true error. In Fig. 7.29, only the error estimators using two or fewer mesh
levels can be used, as the coarsest mesh that could be computed was 33, below which the
solution was found to be numerically unstable. Here both MNP and Richardson
extrapolation using the global order of accuracy gives a slightly higher estimate (but in
reasonable limit). Due to the same reason, in Fig 7.30, only MNP could be used as all
other estimators needed additional coarse meshes. MNP gives a slightly higher estimate
of the error.
x
D
i
s
c
r
et
i
z
at
i
o
n
E
r
r
o
r
-2 -1 0 1 2
-0.0003
-0.0002
-0.0001
0
0.0001
0.0002
0.0003
True Error
MNP
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.27 Discretization error estimators for Burgers equation, Re=64, finest mesh =1025
59
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-2 -1 0 1 2
-0.005
-0.004
-0.003
-0.002
-0.001
0
0.001
0.002
0.003
0.004
0.005
True Error
MNP
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.28 Discretization error estimators for Burgers equation, Re=64, finest mesh =257
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-2 -1 0 1 2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
True Error
MNP
RDE withglobalp
Fig 7.29 Discretization error estimators for Burgers equation, Re=64, finest mesh =65
60
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-1
-0.75
-0.5
-0.25
0
0.25
0.5
0.75
1
True Error
MNP
Fig 7.30 Discretization error estimators for Burgers equation, Re=64, finest mesh =33
We now turn our attention to the discretization error in the nearby problem. Again
different meshes are used for each Reynolds number case. For a Reynolds number of 8,
we used 17 spline fit points. Fig 7.31 shows discretization error estimates with the finest
mesh being 1025 points. In this case, we see that all the error estimators perform well. In
Fig. 7.32, the finest mesh used was 257 points. Again all the estimators give a good
estimate of the discretization error. In Fig. 7.33, the finest mesh that was used was 129
points. Here we see that the mixed-order error estimator begins to under estimate the
error. In Fig. 7.34, both the mixed-order error estimator and Richardson extrapolation
using the local order of accuracy underestimate the error.
61
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-3E-06
-2E-06
-1E-06
0
1E-06
2E-06
3E-06
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.31 Discretization error estimators for the nearby problem, Re=8, finest mesh=1025
x
D
i
s
c
r
et
i
z
at
i
o
n
E
r
r
o
r
-4 -2 0 2 4
-5E-05
-4E-05
-3E-05
-2E-05
-1E-05
0
1E-05
2E-05
3E-05
4E-05
5E-05
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.32 Discretization error estimators for the nearby problem, Re=8, finest mesh=257
62
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-0.0002
-0.00015
-0.0001
-5E-05
0
5E-05
0.0001
0.00015
0.0002
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.33 Discretization error estimators for the nearby problem, Re=8, finest mesh=129
x
D
i
s
c
r
et
i
z
at
i
o
n
E
r
r
o
r
-4 -2 0 2 4
-0.0008
-0.0006
-0.0004
-0.0002
0
0.0002
0.0004
0.0006
0.0008
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.34 Discretization error estimators for the nearby problem, Re=8, finest mesh=65
63
For the Reynolds number of 64 case, we used 33 spline points. Fig. 7.35 shows
the behavior of the error estimators with the finest mesh being 1025 points. In this case,
we see that all the error estimators perform well. In Fig. 7.36, the finest mesh used was
257 points. The mixed-order error estimator and Richardson extrapolation with the local
order of accuracy underestimate the error. In Fig. 7.37, the finest mesh that was used was
129 points. Here the mixed-order error estimator gives very poor estimates and
Richardson extrapolation with the local order of accuracy also begins to fail. In Fig. 7.38,
only Richardson extrapolation using the global order is used, because the other methods
need three mesh solutions, but we cannot go coarser as the numerical solutions become
unstable.
x
D
i
s
c
r
et
i
z
at
i
o
n
E
r
r
o
r
-2 -1 0 1 2
-0.0003
-0.0002
-0.0001
0
0.0001
0.0002
0.0003
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.35 Discretization error estimators for the nearby problem, Re=64, finest
mesh=1025
64
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-2 -1 0 1 2
-0.005
-0.004
-0.003
-0.002
-0.001
0
0.001
0.002
0.003
0.004
0.005
True Error
MixedOrder
RDE withlocalp
RDE withglobalp
Fig 7.36 Discretization error estimators for the nearby problem, Re=64, finest mesh=257
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-0.04
-0.035
-0.03
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.37 Discretization error estimators for the nearby problem, Re=64, finest mesh=129
65
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-2 -1 0 1 2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
True Error
RDE with global p
Fig 7.38 Discretization error estimators for the nearby problem, Re=64, finest mesh=65
If we compare the performance of different error estimators on original Burgers
equation and the nearby problem to Burgers equation, we see that they are very similar.
For example, for Re=8 with finest mesh being 65 points, the error estimator on the
original Burgers equation (Fig. 7.25) are quite similar to the estimates on the nearby
problem (Fig. 7.34). This is also true for the Re=64 case for original Burgers equation
(Fig. 7.28) and the nearby problem (Fig. 7.36). This similarity justifies evaluating error
estimators on the nearby problem.
The nearby problem to the modified form of the Burgers equation was solved and
different error estimators were used to estimate the discretization error. In the modified
form of Burgers equation, we ran only the Re = 64 case. Examining at Figs. 7.39 and
7.40, we see that as we reduce the number of nodes, Richardson extrapolation using local
order of accuracy and the mixed order error estimator breaks down, while Richardson
extrapolation using the global order of accuracy gives the best estimate. In Fig. 7.40,
66
mixed-order error estimator fails while Richardson extrapolation using the local order
also gives bad estimates.
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
-0.001
0
0.001
0.002
0.003
0.004
0.005
0.006
0.007
0.008
0.009
0.01
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.39 Discretization error estimators for nearby problem to modified Burgers equation,
Re=64, nodes=257
x
D
i
s
c
r
e
t
i
z
a
ti
on
E
r
r
o
r
-4 -2 0 2 4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
True Error
MixedOrder
RDE withlocalp
RDE with global p
Fig 7.40 Discretization error estimators for nearby problem to modified Burgers equation,
Re=64, nodes=129
67
CHAPTER EIGHT
CONCLUSIONS AND RECOMMENDATIONS
8.1. Conclusions
The current work on MNP has focused on two major goals. Using MNP to
generate exact solutions to a problem nearby the steady-state Burgers equation in order to
evaluate error estimators was the first goal. The second goal was to use MNP itself as an
error estimator. To fulfill the first goal, a fifth-order Hermite spline was shown to be the
best available curve fitting tool for Burgers equation. The polynomials fits were
demonstrated to be accurate, and we were thus able to generate very small source terms.
MNP was then used to evaluate various discretization error estimators and we found that
Richardson extrapolation using the global order of accuracy provided the best estimates.
MNP itself was used as an error estimator. It was shown that MNP provided error
estimates that were at least as good as Richardson extrapolation with the global order of
accuracy, and often better. The use of MNP as an error estimator requires numerical
solutions on only one mesh. This is an advantage of using MNP, as other error estimators
required multiple meshes. The cost of MNP alone will be approximately the same as the
cost of solving the original problem, with some additional overhead involved in the curve
fitting procedure. So the total cost of using MNP as an error estimator will be
approximately twice the cost of numerically solving the original problem of interest
alone.
68
8.2. Future work
This research was focused on one-dimensional problems. The next step is to
extend the MNP methodology to two-dimensions, three-dimensions, and four-dimensions
(three spatial dimensions plus time). The biggest challenge is to develop a two-
dimensional or higher-dimensional fit which has continuity. Rouff [17] has developed
an approach for generating n-dimensional, continuous splines. Implementation of
Rouff?s approach will allow the extension of MNP to more realistic two-dimensional and
three-dimensional applications.
3
C
k
C
69
REFERENCES
1. Anderson, J. D., ?Computational Fluid Dynamics: The Basics With Applications,?
McGraw-Hill, New York, 1995
2. Oberkampf, W. L., Roy, C. J., short course for Moldflow corporation on ?Verification
and Validation in Computational Simulation,? Jul 22-23, 2004, Ithaca, New York
3. Roy, C. J., ?Review of Code and Solution Verification Procedures for Computational
Simulation? Journal of Computational Physics, Vol. 205, 2005, pp. 131-156
4. www.softwareqatest.com, access date: June 29, 2005
5. Ferziger, J. H., Peric, M., ?Computational Methods of Fluid Dynamics,? Springer-
Verlag, Berlin, 1999
6. Raju, A., Roy, C. J., Hopkins, M. M., ?Evaluation of Discretization Error Estimators
using the Method of Nearby Problems,? AIAA 2005-0684, 35
th
AIAA Fluid
Dynamics Conference & Exhibit, Toronto, Canada
7. Roache, P. J., ?Code Verification by the Method of Manufactured Solutions,? Journal
of Fluids Engineering. 124(1) (2002) 4-10
8. Knupp, P., Salari, K., ?Verification of Computer Codes in Computational Science and
Engineering,? Chapman & Hall/CRC, New York, 2002
9. Lee, S., Junkins, J. L., ?Construction of Benchmark Problems for Solution of ordinary
Differential Equations,? Journal of Shock and Vibration, Vol. 1, No. 5, 1994, pp.
403-414
70
10. Junkins, J. L, Lee, S., ?Validation of Finite-Dimensional Approximate Solutions for
Dynamics of Distributed-Parameter Systems,? Journal of Guidance, Control, and
Dynamics, Vol. 18, No. 1, 1995, pp. 87-95
11. Roy, C. J., and Hopkins, M. M., ?Discretization Error Estimates using Exact
Solutions to Nearby Problems,? AIAA Paper 2003-0629, January 2003
12. Benton, E. R., and Platzman, G. W., ?A Table of Solutions of the One-Dimensional
Burger?s Equation,? Q. Appl Math, Vol. 30. pp. 195-212, 1972
13. Hopkins, M. M., Roy, C. J., ?Introducing The Method of Nearby Problems,?
European Congress on Computational Methods in Applied Sciences and Engineering,
ECCOMAS 2004, P. Neittaanmaki, T. Rossi, S. Korotov, E. Onate, J. Periaux, and D.
Knorzer (eds.), July 24-28, 2004
14. Matlab user manual
15. Mullges, G. E., Uhlig, F., ?Numerical Algorithms with Fortran,? Springer-Verlag,
Berlin, 1996
16. Roy, C. J., ?Grid Convergence Error Analysis for Mixed-Order Numerical Schemes,?
AIAA journal, Vol. 41, No. 4, April 2003
17. Rouff, M., ?The Computation of Ck Spline Functions,? Computers in Mathematical
Applications, Vol 23, No. 1, 1992, pp. 103-110
71
APPENDICES
A. Fortran program for solving original Burgers equation
Program Steady_implicit
implicit doubleprecision(a-h,o-z)
Parameter(imax=81,lmax=imax-2)
Dimension uold(imax),up(imax),unew(imax),x(imax)
Dimension
uexact(imax),Disc_error(imax),AA(imax),BB(imax),CC(imax)
Dimension G(imax),y(imax),ub(imax),unew1(imax),unew2(imax)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxx Initialization of constants xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
rLen = 4.0d0
a=-4.0d0
b=4.0d0
rNu=2.0d0
rH=rLen*2.0d0/dfloat(imax-1)
cfl=200.0d0
rK=cfl*((rH*rH)/(rH + 2.0d0*rNu))
n=imax-1
th=1.0d0
t=0.0d0
Tol = 2.7e-12
jmax = 400000000
al=2.0d0
alp=rK/(2.0d0*rH)
bet=rNu*rK/rH**2
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Initial Conditions xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i = 1,imax
x(i) = a + (b-a)*dfloat(i-1)/dfloat(n)
! uold(i)=0.0
uold(i) = (-2.0d0*rNu*al/rLen)*(dsinh(x(i)*al/rLen))/
& (dcosh(x(i)*al/rLen)+dexp(-(t*rNu*al**2/rLen**2)))
! Write(*,*) x(i)
End Do
72
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Boundary Conditions xxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
uold(1) = -2.0d0*al*rNu/rLen*dtanh(a*al/rLen)
uold(n+1)=-2.0d0*al*rNu/rLen*dtanh(b*al/rLen)
Do i=1,imax
uexact(i)=-2.0d0*rNu*al/rLen*dtanh(x(i)*al/rLen)
End Do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxxxx Main Loop xxxxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do j=1,jmax
Do i=1,imax
ub(i)=uold(i)
End do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Setting up of Tridiagonal Matrix xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=3,imax-2
k=i-1
AA(k)=-1.0d0*th*(alp*ub(i)+bet)
BB(k)=(1.0d0+2.0d0*th*bet)
CC(k)=th*(alp*ub(i)-bet)
G(k)=uold(i)-((1.0d0-th)*alp*ub(i)*(uold(i+1)-uold(i-1)))+
& ((1.0d0-th)*bet*(uold(i+1)-2.0d0*uold(i)+uold(i-1)) )
End Do
BB(1)=(1.0d0+2.0d0*th*bet)
CC(1)= (th*alp*ub(2)-th*bet)
AA(imax-2)=-1.0d0*(th*alp*ub(imax-1)+th*bet)
BB(imax-2)=(1.0d0+2.0d0*th*bet)
G(1)=uold(2)-((1.0d0-th)*alp*ub(2)*(uold(3)-uold(1)))+
& ((1.0d0-th)*bet*(uold(3)-2.0d0*uold(2)+uold(1)))+(th*alp*
& ub(2)+th*bet)*uold(1)
G(imax-2)=uold(imax-1)-((1.0d0-th)*alp*ub(imax-1)*(uold(imax)-
& uold(imax-2)))+((1.0d0-th)*bet*(uold(imax)-2.0d0*uold(imax-
1)+
& uold(imax-2)))-(th*alp*ub(imax-1)-th*bet)*uold(imax)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Calling the Tridiagonal Solver xxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
CALL TRDIAG (lmax,AA,BB,CC,y,G)
73
Do i=1,imax-2
unew(i+1)=y(i)
End do
unew(1)=uold(1)
unew(n+1)=uold(n+1)
Residue = 0.0d0
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxx Computation of Residual xxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=2,imax-1
Res = unew(i)/(2.0d0*rH)*(unew(i+1)-unew(i-1))-rNu/(rH**2)*
& (unew(i+1)-2.0d0*unew(i)+unew(i-1))
Residue = Residue + Res**2
End Do
rNorm=Residue
rL2Norm =dsqrt(rNorm/dfloat(imax-2))
if (j.eq.1) r1 = rL2Norm
rL2Norm=rL2Norm/r1
write(11,*) j,rK*dfloat(j),rK,rL2Norm
if(mod(j,10).eq.0.0d0) then
Write(19,*) j,rL2Norm
End if
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxx Checking Convergence xxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
If (rL2Norm.LT.Tol) Goto 100
200 Do i =1,imax
uold(i) = unew(i)
End Do
End Do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxx End of Main Loop xxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
100 Write(*,*) j
do i=2,imax-1
unew1(i)=(unew(i+1)-unew(i-1))/(2.0d0*rH)
End do
unew1(1)=(-3.0d0*unew(1)+4.0d0*unew(2)-unew(3))/(2.0d0*rH)
unew1(imax)=(3.0d0*unew(imax)-4.0d0*unew(imax-1)+unew(imax-2))/
& (2.0d0*rH)
unew2(1)=(-unew(4)+4.0d0*unew(3)-5.0d0*unew(2)+2.0d0*unew(1))/
& rH**2
unew2(imax)=(2.0d0*unew(imax)-5.0d0*unew(imax-1)+4.0d0*
74
& unew(imax-2)-unew(imax-3))/rH**2
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Writing out Solution xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
open(30,file='velocity_der_new_17.dat',status='unknown')
open(40,file='distance_new_17.dat',status='unknown')
open(50,file='velocity_new_17.dat',status='unknown')
Do i=1, imax,5
write(30,6) unew1(i)
write(40,6) x(i)
write(50,6) unew(i)
6 Format (e30.20)
end do
open(15,file='velocity.dat',status='unknown')
open(10,file='distance1.dat',status='unknown')
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxx Computation of Discretization Error xxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=1,imax
Write (10,500) x(i)
Write (15,500) unew(i)
500 Format(e30.20)
Write (16,501) x(i),unew(i),uexact(i)
Disc_error(i)=(unew(i)-uexact(i))
Write(14,*) x(i),Disc_error(i)
501 Format(3(e30.20))
Write(99,500) uexact(i)
End Do
rl2 = 0.0d0
rl1 = 0.0d0
Do i=1,imax
error = dabs(uexact(i) - unew(i))
error2 = error*error
rl1 = rl1 + error
rl2 = rl2 + error2
End Do
rl1 = rl1/dfloat(n+1)
rl2 = dsqrt(rl2/dfloat(n+1))
Write(*,*) 'L2Norm=',rl2
Write(*,*) 'L1Norm=',rl1
Write(*,8) unew2(1)
Write(*,8) unew2(imax)
8 Format(e30.20)
close(30)
close(40)
close(50)
End
75
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Tridiagonal Solver Subroutine xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
SUBROUTINE TRDIAG (N,A,B,C,X,G)
implicit doubleprecision(a-h,o-z)
DIMENSION A(2000),B(2000),C(2000),X(2000),G(2000),BB(2000)
C!.....THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF EQUATIONS
C!.....BY GAUSS ELIMINATION
C!.....THE PROBLEM SOLVED IS MX=G WHERE M=TRI(A,B,C)
C!.....THIS ROUTINE DOES NOT DESTROY THE ORIGINAL MATRIX
C!.....AND MAY BE CALLED A NUMBER OF TIMES WITHOUT REDEFINING
C!.....THE MATRIX
C!.....N = NUMBER OF EQUATIONS SOLVED (UP TO 1000)
C!.....FORWARD ELIMINATION
C!.....BB IS A SCRATCH ARRAY NEEDED TO AVOID DESTROYING B ARRAY
DO 1 I=1,N
BB(I) = B(I)
1 CONTINUE
DO 2 I=2,N
T = A(I)/BB(I-1)
BB(I) = BB(I) - C(I-1)*T
G(I) = G(I) - G(I-1)*T
2 CONTINUE
C!.....BACK SUBSTITUTION
X(N) = G(N)/BB(N)
DO 3 I=1,N-1
J = N-I
X(J) = (G(J)-C(J)*X(J+1))/BB(J)
3 CONTINUE
RETURN
END
B. Fortran program for solving the nearby problem to Burgers equation
Program Steady_implicit
implicit doubleprecision(a-h,o-z)
Parameter(imax=81)
Dimension uold(imax),unew(imax),x(imax),u(imax)
Dimension
source(imax),Disc_error(imax),AA(imax),BB(imax),CC(imax)
Dimension G(imax),ub(imax),y(imax)
rLen = 4.0d0
10 Format (e25.14)
76
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxx Initialization of constants xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
a=-4.0d0
b=4.0d0
rNu=2.0d0
rH=rLen*2.0d0/dfloat(imax-1)
cfl=100.0d0
rK=cfl*((rH*rH)/(rH + 2.0d0*rNu))
n=imax-1
th=1.0d0
! write(*,*) 'n = ',n
t=0.0d0
Tol = 1e-12
jmax = 400000000
al=2.0d0!255.828300733!7.99463621669!4.0!255.828300733
alp=rK/(2.0d0*rH)
bet=rNu*rK/rH**2
open(60,file='distance1.dat',status='unknown')
open(70,FILE='velocityterm.dat',status='unknown')
open(66,file='source_17.dat',status='unknown')
Do i=1,imax
Read(60,*) x(i)
Read(70,*) u(i)
Read(66,*) source(i)
End do
close(60)
close(70)
close(66)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Initial Conditions xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i = 1,imax
x(i) = a + (b-a)*dfloat(i-1)/dfloat(n)
uold(i) = (-2.0d0*rNu*al/rLen)*(dsinh(x(i)*al/rLen))/
& (dcosh(x(i)*al/rLen)+dexp(-(t*rNu*al**2/rLen**2)))
End Do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Boundary Conditions xxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
uold(1) = u(1)!-2.0d0*al*rNu/rLen*dtanh(a*al/rLen)
uold(n+1)=u(imax)!-2.0d0*al*rNu/rLen*dtanh(b*al/rLen)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxxxx Main Loop xxxxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
77
Do j=1,jmax
Do i=1,imax
ub(i)=uold(i)
End do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Setting up of Tridiagonal Matrix xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=2,imax-1
AA(i)=-1.0d0*th*(alp*ub(i)+bet)
BB(i)=(1.0d0+2.0d0*th*bet)
CC(i)=th*(alp*ub(i)-bet)
G(i)=source(i)*rK+uold(i)-((1.0d0-
th)*alp*ub(i)*(uold(i+1)-
& uold(i-1)))+((1.0d0-th)*bet*(uold(i+1)-2.0d0*uold(i)+
& uold(i-1)) )
End Do
BB(1)=1.0d0
CC(1)= 0.0d0
AA(imax)=0.0d0
BB(imax)=1.0d0
G(1)=u(1)
G(imax)=u(imax)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Calling the Tridiagonal Solver xxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
CALL TRDIAG (imax,AA,BB,CC,y,G)
Do i=1,imax
unew(i)=y(i)
End do
Residue = 0.0d0
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxx Computation of Residual xxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=2,imax-1
Res =unew(i)-uold(i)! unew(i)/(2.0d0*rH)*(unew(i+1)-unew(i-
1))-rNu/(rH**2)*
Residue = Residue + Res**2
End Do
rNorm=Residue
rL2Norm =dsqrt(rNorm/dfloat(imax-2))
78
if (j.eq.1) r1 = rL2Norm
rL2Norm=rL2Norm/r1
write(11,*) j,rK*dfloat(j),rK,rL2Norm
if(mod(j,10).eq.0.0d0) then
Write(19,*) j,rL2Norm
End if
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxx Checking Convergence xxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
If (rL2Norm.LT.Tol) Goto 100
Do i =1,imax
uold(i) = unew(i)
End Do
End Do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxx End of Main Loop xxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
100 Write(*,*) j
open(14,file='MNP.dat',status='unknown')
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxx Computation of Discretization Error xxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=1,imax
Disc_error(i)=(unew(i)-u(i))/u(i)
Write(14,8) x(i),Disc_error(i)
Write (16,9) x(i),unew(i),u(i)
End Do
8 Format(2(e30.20))
9 Format(3(e30.20))
rl2 = 0.0d0
rl1 = 0.0d0
Do i=1,imax
error = dabs(u(i) - unew(i))
error2 = error*error
rl1 = rl1 + error
rl2 = rl2 + error2
End Do
rl1 = rl1/dfloat(imax)
rl2 = dsqrt(rl2/dfloat(imax))
Write(*,99) rl1
79
Write(*,99) rl2
99 Format(e30.20)
End
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Tridiagonal Solver Subroutine xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
SUBROUTINE TRDIAG (N,A,B,C,X,G)
implicit doubleprecision(a-h,o-z)
DIMENSION A(2000),B(2000),C(2000),X(2000),G(2000),BB(2000)
C!.....THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF EQUATIONS
C!.....BY GAUSS ELIMINATION
C!.....THE PROBLEM SOLVED IS MX=G WHERE M=TRI(A,B,C)
C!.....THIS ROUTINE DOES NOT DESTROY THE ORIGINAL MATRIX
C!.....AND MAY BE CALLED A NUMBER OF TIMES WITHOUT REDEFINING
C!.....THE MATRIX
C!.....N = NUMBER OF EQUATIONS SOLVED (UP TO 1000)
C!.....FORWARD ELIMINATION
C!.....BB IS A SCRATCH ARRAY NEEDED TO AVOID DESTROYING B ARRAY
DO 1 I=1,N
BB(I) = B(I)
1 CONTINUE
DO 2 I=2,N
T = A(I)/BB(I-1)
BB(I) = BB(I) - C(I-1)*T
G(I) = G(I) - G(I-1)*T
2 CONTINUE
C!.....BACK SUBSTITUTION
X(N) = G(N)/BB(N)
DO 3 I=1,N-1
J = N-I
X(J) = (G(J)-C(J)*X(J+1))/BB(J)
3 CONTINUE
RETURN
END
C. Fortran program to numerically solve the modified form of Burgers equation
Program Steady_implicit
implicit doubleprecision(a-h,o-z)
Parameter(imax=1025,lmax=imax-2)
Dimension uold(imax),up(imax),unew(imax),x(imax)
Dimension
uexact(imax),Disc_error(imax),AA(imax),BB(imax),CC(imax)
Dimension G(imax),y(imax),ub(imax),unew1(imax),unew2(imax)
80
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxx Initialization of constants xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
rLen = 4.0d0
a=-4.0d0
b=4.0d0
rNu=0.25d0 !2.0d0, 0.25d0, 0.03125d0
rH=rLen*2.0d0/dfloat(imax-1)
cfl=1000.0d0
rK=cfl*((rH*rH)/(rH + 2.0d0*rNu))
n=imax-1
th=1.0d0
t=0.0d0
Tol = 3e-12
jmax = 400000000
al=16.0d0
alp=rK/(2.0d0*rH)
bet=rNu*rK/rH**2
xxp = 0.25
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Boundary Conditions xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
x(1) = a
x(imax) = b
uold(1) = 2.
uold(n+1) = -2.
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Setting up of Initial Conditions xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i = 2,imax-1
x(i) = a + (b-a)*dfloat(i-1)/dfloat(n)
uold(i) = uold(1)+(uold(n+1)-uold(1))*dfloat(i-1)/dfloat(n)
End Do
Write(11,*) 'TITLE = "Solution Profiles"'
Write(11,*) 'variables="Iterations""rk*float(j)""rK""L2Norm"'
Do i=1,imax
uexact(i)=-2.0d0*rNu*al/rLen*dtanh(x(i)*al/rLen)
End Do
81
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxxxx Main Loop xxxxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do j=1,jmax
Do i=1,imax
ub(i)=uold(i)
End do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Setting up of Tridiagonal Matrix xxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=3,imax-2
k=i-1
betvar=bet*( (ub(i)/2.)**2 + ( (x(i)-a)/(b-a)+0.25)**xxp )
AA(k)=-1.0d0*th*(alp*ub(i)+betvar)
BB(k)=(1.0d0+2.0d0*th*betvar)
CC(k)=th*(alp*ub(i)-betvar)
G(k)=uold(i)-((1.0d0-th)*alp*ub(i)*(uold(i+1)-uold(i-1)))+
& ((1.0d0-th)*betvar*(uold(i+1)-2.0d0*uold(i)+uold(i-1)) )
End Do
betvar=bet*( (ub(2)/2.)**2 + ( (x(2)-a)/(b-a) + 0.25)**xxp )
BB(1)=(1.0d0+2.0d0*th*betvar)
CC(1)= (th*alp*ub(2)-th*betvar)
G(1)=uold(2)-((1.0d0-th)*alp*ub(2)*(uold(3)-uold(1)))+
& ((1.0d0-th)*betvar*(uold(3)-2.0d0*uold(2)+uold(1)))+(th*alp*
& ub(2)+th*betvar)*uold(1)
betvar=bet*( (ub(imax-1)/2.)**2 +
& ( (x(imax-1)-a)/(b-a) + 0.25)**xxp )
AA(imax-2)=-1.0d0*(th*alp*ub(imax-1)+th*betvar)
BB(imax-2)=(1.0d0+2.0d0*th*betvar)
G(imax-2)=uold(imax-1)-
& ((1.0d0-th)*alp*ub(imax-1)*(uold(imax)-
& uold(imax-2)))+((1.0d0-th)*betvar*(uold(imax)
& -2.0d0*uold(imax-1)+
& uold(imax-2)))-(th*alp*ub(imax-1)-th*betvar)*uold(imax)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Calling the Tridiagonal Solver xxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
CALL TRDIAG (lmax,AA,BB,CC,y,G)
Do i=1,imax-2
unew(i+1)=y(i)
End do
unew(1)=uold(1)
unew(n+1)=uold(n+1)
Residue = 0.0d0
82
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxx Computation of Residual xxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=2,imax-1
rnuvar=rNu*( (unew(i)/2.)**2 + ( (x(i)-a)/(b-a)+0.25)**xxp )
Res = unew(i)/(2.0d0*rH)*(unew(i+1)-unew(i-1))-
& rnuvar/(rH**2)*(unew(i+1)-2.0d0*unew(i)+unew(i-1))
Residue = Residue + Res**2
End Do
rNorm=Residue
rL2Norm =dsqrt(rNorm/dfloat(imax-2))
if (j.eq.1) r1 = rL2Norm
rL2Norm=rL2Norm/r1
write(11,*) j,rK*dfloat(j),rK,rL2Norm
if(mod(j,100).eq.0.0d0) then
Write(19,*) j,rL2Norm
write(*,*) j,rL2Norm
End if
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxx Checking Convergence xxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
If (rL2Norm.LT.Tol) Goto 100
200 Do i =1,imax
uold(i) = unew(i)
End Do
End Do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxxxxx End of Main Loop xxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
100 Write(*,*) j
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Writing out Solution xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
do i=2,imax-1
unew1(i)=(unew(i+1)-unew(i-1))/(2.0d0*rH)
End do
unew1(1)=(-3.0d0*unew(1)+4.0d0*unew(2)-unew(3))/(2.0d0*rH)
unew1(imax)=(3.0d0*unew(imax)-4.0d0*unew(imax-1)+unew(imax-2))/
& (2.0d0*rH)
unew2(1)=(-unew(4)+4.0d0*unew(3)-5.0d0*unew(2)+2.0d0*unew(1))/
& rH**2
unew2(imax)=(2.0d0*unew(imax)-5.0d0*unew(imax-1)+4.0d0*
& unew(imax-2)-unew(imax-3))/rH**2
open(30,file='velocity_der_new_33_F.dat',status='unknown')
open(40,file='distance_new_33_F.dat',status='unknown')
83
open(50,file='velocity_new_33_F.dat',status='unknown')
Do i=1, imax,32
write(30,6) unew1(i)
write(40,6) x(i)
write(50,6) unew(i)
6 Format (e30.20)
end do
open(15,file='velocity_F.dat',status='unknown')
open(10,file='distance1_F.dat',status='unknown')
open(51,file='distance1_M.dat',status='unknown')
open(52,file='distance1_C.dat',status='unknown')
open(11,file='ubar_F.dat',status='unknown')
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxx Computation of Discretization Error xxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=1,imax
Write (10,500) x(i)
Write (15,500) unew(i)
write(11,500) ub(i)
500 Format(e30.20)
Disc_error(i)=(unew(i)-uexact(i))
Write(14,*) x(i),Disc_error(i)
End Do
Do i=1,imax,2
Write (51,501) x(i)
501 Format(e30.20)
Enddo
Do i=1,imax,4
Write (52,502) x(i)
502 Format(e30.20)
Enddo
rl2 = 0.0d0
rl1 = 0.0d0
Do i=1,imax
error = dabs(uexact(i) - unew(i))
error2 = error*error
rl1 = rl1 + error
rl2 = rl2 + error2
End Do
rl1 = rl1/dfloat(n+1)
rl2 = dsqrt(rl2/dfloat(n+1))
write(*,*) rl1,rl2
84
Write(*,8) rl2
! Write(*,*) 'L1Norm=',rl1
Write(*,8) unew2(1)
Write(*,8) unew2(imax)
8 Format(e30.20)
close(30)
close(40)
close(50)
End
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxx Tridiagonal Solver Subroutine xxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
SUBROUTINE TRDIAG (N,A,B,C,X,G)
implicit doubleprecision(a-h,o-z)
DIMENSION A(2000),B(2000),C(2000),X(2000),G(2000),BB(2000)
C!.....THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF EQUATIONS
C!.....BY GAUSS ELIMINATION
C!.....THE PROBLEM SOLVED IS MX=G WHERE M=TRI(A,B,C)
C!.....THIS ROUTINE DOES NOT DESTROY THE ORIGINAL MATRIX
C!.....AND MAY BE CALLED A NUMBER OF TIMES WITHOUT REDEFINING
C!.....THE MATRIX
C!.....N = NUMBER OF EQUATIONS SOLVED (UP TO 1000)
C!.....FORWARD ELIMINATION
C!.....BB IS A SCRATCH ARRAY NEEDED TO AVOID DESTROYING B ARRAY
DO 1 I=1,N
BB(I) = B(I)
1 CONTINUE
DO 2 I=2,N
T = A(I)/BB(I-1)
BB(I) = BB(I) - C(I-1)*T
G(I) = G(I) - G(I-1)*T
2 CONTINUE
C!.....BACK SUBSTITUTION
X(N) = G(N)/BB(N)
DO 3 I=1,N-1
J = N-I
X(J) = (G(J)-C(J)*X(J+1))/BB(J)
3 CONTINUE
RETURN
END
85
D. Fortran program to compute the coefficients of the Hermite spline polynomial
Program Spline
implicit doubleprecision(a-h,o-z)
Parameter (imax=257,kmax=2000,nmax=33,lmax=nmax-2)
Dimension x(kmax),u(kmax),a(kmax),b(kmax),c(kmax),d(kmax),h(kmax)
Dimension
g(kmax),aa(kmax),bb(kmax),cc(kmax),u_num(kmax),u_1(kmax)
Dimension e(kmax),f(kmax),y(kmax),z(kmax),p(kmax),q(kmax)
Dimension term(kmax),de2(imax)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Reading Input Files xxxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
open(60,file='distance_new_33.dat',status='unknown')
open(70,FILE='velocity_new_33.dat',status='unknown')
open(80,FILE='velocity_der_new_33.dat',status='unknown')
open(66,FILE='distance1.dat',status='unknown')
open(55,FILE='velocity.dat',status='unknown')
Do i=1,nmax
Read(60,*) x(i)
Read(70,*) u(i)
Read(80,*) u_1(i)
!100 Format (F7.4)
End do
Do i=1,imax
Read(55,*) u_num(i)
Read(66,*) y(i)
Enddo
close(60)
close(70)
close(80)
close(66)
close(55)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxx Initializing Constants xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
toler=0.25d0
rh=0.25d0
rnu=0.250d0
alp=1.0d0
bet1=-0.5d0*0.90949470177292824000E-12/rh
bet2=-0.5d0*0.20463630789890885000E-11/rh
86
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxx Initializing ?a? and ?b? Coefficients xxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
do i = 1,nmax
a(i)=u(i)
b(i)=u_1(i)
end do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxx Setting up Tridiagonal Matrix for the Coefficient ?c? xxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
c(1)=-0.5d0*0.90949470177292824000E-12!0.5/h(1)*(3.0/h(1)*(a(2)-
a(1))-3.0*ral1-ral1/2.0
c(nmax)=-0.5d0*0.20463630789890885000E-11!0.5*0.0101!ral2/2.0
i=2
k=i-1
g(k)=10.0d0/rh**3*(a(3)-2.0d0*a(2)+a(1))+4.0d0/rh**2*(b(1)-b(3))+
& 1.0d0/rh*c(1)
i=nmax-1
k=i-1
g(k)= 10.0d0/rh**3*(a(nmax)-2.0d0*a(nmax-1)+a(nmax-2))+
& 4.0d0/rh**2*(b(nmax-2)-b(nmax))+1.0d0/rh*c(nmax)
do i=3,nmax-2
k=i-1
g(k)=10.0d0/rh**3*(a(i+1)-2.0d0*a(i)+a(i-1))+4.0d0/rh**2*
& (b(i-1)-b(i+1))
aa(k)= -1.0d0/rh
bb(k)= 6.0d0/rh
cc(k)= -1.0d0/rh
end do
bb(1)=6.0d0/rh
cc(1)=-1.0d0/rh
aa(lmax)=-1.0d0/rh
bb(lmax)=6.0d0/rh
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxx Calling Tridiagonal Solver xxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
CALL TRDIAG (lmax,aa,bb,cc,z,g)
Do i=1,nmax-2
87
c(i+1)=z(i)
End do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxx Computing?d?, ?e?, and ?f? Coefficients xxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=1,nmax-1
d(i)=10.0d0/rh**3*(a(i+1)-a(i))-2.0d0/rh**2*(2.0d0*b(i+1)+
& 3.0d0*b(i))+1.0d0/rh*(c(i+1)-3.0d0*c(i))
e(i)=5.0d0/rh**4*(a(i+1)-a(i))-
1.0d0/rh**3*(b(i+1)+4.0d0*b(i))
& -3.0d0/rh**2*c(i)-2.0d0/rh*d(i)
f(i)=1.0d0/(10.0d0*rh**3)*(c(i+1)-c(i)-3.0d0*d(i)*rh-6.0d0*
& e(i)* rh**2)
End do
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxx Writing out Coefficients xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
open(20,file='OUTPUT_hermite_new_33.dat',status='unknown')
Do i=1,nmax-1
write(20,200) a(i), b(i), c(i), d(i), e(i), f(i)
200 Format(6(e30.20))
End do
close(20)
open(77,file='G.dat',status='unknown')
open(88,file='H.dat',status='unknown')
Do i=1,imax-1
Read(77,*) p(i)
Read(88,*) q(i)
Write(*,*) p(i),q(i)
End do
close(77)
close(88)
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxx Computing Value of the Polynomial xxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
Do i=1,imax
Do j=1,nmax-1
if((y(i)-x(j)).le.toler) goto 300
enddo
300 term(i)=a(j)+b(j)*(y(i)-x(j))+c(j)*(y(i)-x(j))**2+d(j)*
& (y(i)-x(j))**3+e(j)*(y(i)-x(j))**4+f(j)*(y(i)-x(j))**5
enddo
88
rl1norm=0.0d0
Do i=1,imax-1
Do j=1,nmax
if ((y(i)-x(j)).le.toler) goto 400
enddo
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Computing the Norm xxxxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
400 rl1norm = rl1norm + dabs(((p(i)*y(i+1)+q(i)/2.0d0*y(i+1)**2)-
& (a(j)*y(i+1)+b(j)/2.0d0*(y(i+1)-x(j))**2+c(j)/3.0d0*(y(i+1)-
& x(j))**3+d(j)/4.0d0*(y(i+1)-x(j))**4+e(j)/5.0d0*(y(i+1)-
& x(j))**5+f(j)/6.0d0*(y(i+1)-x(j))**6))-((p(i)*y(i)+q(i)/
& 2.0d0*y(i)**2)-(a(j)*y(i)+b(j)/2.0d0*(y(i)-x(j))**2+c(j)/
& 3.0d0*(y(i)-x(j))**3+d(j)/4.0d0*(y(i)-x(j))**4+e(j)/5.0d0*
& (y(i)-x(j))**5+f(j)/6.0d0*(y(i)-x(j))**6)))
end do
rl1norm=rl1norm/8
Write(*,145) rl1norm
145 Format(e30.20)
open(22,file='term_33.dat',status='unknown')
open(44,file='velocityterm.dat',status='unknown')
Do i=1,imax
Write(44,141) term(i)
141 format(e30.20)
Write(22,140) y(i),term(i),u_num(i)
140 Format (3(e30.20))
End do
Do i=1,imax
de2(i)=u_num(i)-term(i)
End do
open(23,file='de2.dat',status='unknown')
Do i=1,imax
Write(23,333) y(i),term(i),u_num(i),de2(i)
333 Format (4(e30.20))
End do
close(22)
end
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxx Tridiagonal Solver Subroutine xxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
SUBROUTINE TRDIAG (N,A,B,C,X,G)
implicit doubleprecision(a-h,o-z)
DIMENSION A(2000),B(2000),C(2000),X(2000),G(2000),BB(2000)
C!.....THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF EQUATIONS
89
C!.....BY GAUSS ELIMINATION
C!.....THE PROBLEM SOLVED IS MX=G WHERE M=TRI(A,B,C)
C!.....THIS ROUTINE DOES NOT DESTROY THE ORIGINAL MATRIX
C!.....AND MAY BE CALLED A NUMBER OF TIMES WITHOUT REDEFINING
C!.....THE MATRIX
C!.....N = NUMBER OF EQUATIONS SOLVED (UP TO 1000)
C!.....FORWARD ELIMINATION
C!.....BB IS A SCRATCH ARRAY NEEDED TO AVOID DESTROYING B ARRAY
DO 1 I=1,N
BB(I) = B(I)
1 CONTINUE
DO 2 I=2,N
T = A(I)/BB(I-1)
BB(I) = BB(I) - C(I-1)*T
G(I) = G(I) - G(I-1)*T
2 CONTINUE
C!.....BACK SUBSTITUTION
X(N) = G(N)/BB(N)
DO 3 I=1,N-1
J = N-I
X(J) = (G(J)-C(J)*X(J+1))/BB(J)
3 CONTINUE
RETURN
END
E. Matlab program to calculate the source terms
clc
clear all
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Reading in Coefficients xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
op=load('OUTPUT_hermite_new_17.dat');
y=load('distance1.dat');
x=load('distance_new_17.dat');
a=op(:,1);
b=op(:,2);
c=op(:,3);
d=op(:,4);
e=op(:,5);
f=op(:,6);
tol=0.50;
nu=2.0;
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxx Computation of Source Term xxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
for i=1:81
for j=1:16
if y(i)-x(j)<=tol ,break,end
end
90
source(i)=(a(j)+b(j)*(y(i)-x(j))+c(j)*(y(i)-x(j))^2+d(j)*(y(i)-
x(j))^3+e(j)*(y(i)-x(j))^4+f(j)*(y(i)-x(j))^5)*(b(j)+2*c(j)*(y(i)-
x(j))+3*d(j)*(y(i)-x(j))^2+4*e(j)*(y(i)-x(j))^3+5*f(j)*(y(i)-x(j))^4)-
nu*(2*c(j)+6*d(j)*(y(i)-x(j))+12*e(j)*(y(i)-x(j))^2+20*f(j)*(y(i)-
x(j))^3);
end
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
xxxxxxxxxxxxxxxxxxxxx Plotting Source Term xxxxxxxxxxxxxxxxxxxxxx
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
source=source';
plot(y,source);
l1=0;
l2=0;
for i=1:81
res2=source(i)^2;
res1=abs(source(i));
l1=l1+res1;
l2=l2+res2;
end
l1=l1/81;
l2=(l2/81)^0.5;
91