Teacher Manual
Teacher Manual
Teacher Manual
INTRODUCTION
The M.Sc. Previous course on Computer Programming is meant to introduce students to both com-
puter programming as well as to numerical analysis. Although it expects students not to have forgotten
what they learnt about numerical methods in their B.Sc. Honours course, it makes no great demands
on them in this respect.
The name Linux was coined in 1991 by Linus Torvalds after his own name and is supposed to indicate
any open sourced operating system based on a Unix platform. The word open source implies that
typically all underlying source codes can be freely modified, used, and redistributed by anyone. The
system’s utilities and libraries usually come from the GNU operating system, announced in 1983 by
Richard Stallman.
There are several distributions (popularly referred to as distros) of the Linux OS eg. Red Hat, SUSE
Linux (Novell) and Debian; and other derived distros like Fedora, Ubuntu (Canonical Ltd.), Gentoo,
PCLinuxOS, CentOS, etc. Broadly speaking, the mother distros (Red Hat, SuSE, Debian) differ mostly
in the areas they specialize or otherwise their focus.
In this course we shall use Fedora based on Red Hat Enterprise distribution, Version:
Fedora
1
INTRODUCTION COMPUTER PROGRAMMING-2015
1.1.2 COMPILER
The name“compiler” is primarily used for programs that translate source codes from a high-level pro-
gramming language (like C, C++, Fortran, etc.) to a lower level language (e.g., assembly language or
machine language). The original sequence is usually called the source code and the output the object
code. The most common reason for wanting to translate source codes is to create an executable program.
Some useful linux commands are given here (you dont have to type $, which shows the command
prompt).
Page 2 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
$pwd
print working directory. The command enlists in the present work directory.
$ls
lists files in a given directory. There are many options available with ls, which you can find using the
command man ls or ls –help.
$cd
change directory. To see what is stored in any subdirectory, we must first of all change the current
directory by using the cd command. Suppose we wish to change to the include directory. The required
command then is
$cd include
This when followed by the ls command will display the files and subdirectories stored in the include
directory. To return to the main directory we must execute
$cd ..
Here .. means the parent directory (or one directory up).
$mkdir mydirectory
creates a directory with the name mydirectory
$rm -i filename1.c
deletes/removes a file. Once again the flag -i ensures a check before the command is executed.
$gcc -- version
This will type the version of the compiler you are using. Note that this compiler will compile programs
written in C.
The language C is case sensitive. For example, the variable lower case x is different from
the variable upper case X in C.
Page 3 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
1. First, you need to write and save your program in a text file using a TEXT EDITOR
2. Then you compile this program file using the gcc compiler to create an executable file.
3. Finally, you can EXECUTE or RUN the executable file created in step 2 to get the output
You need a text editor to type and save your program. In Linux there are several editors, like gedit,
emacs, vim etc. You can use any one of these to make your first program. Here, we are showing you
how to make your first program using emacs, but it can be done using any of your favorite editors. On
the command prompt of a terminal, type
$emacs &
This will open a new screen. Using keyboard, type the following program:
Program 1.1
/∗ This program p r i n t s ‘ ‘My F i r s t C programme ’ ’ ∗/
#i n c l u d e <s t d i o . h>
#i n c l u d e <math . h>
main ( )
{
p r i n t f ( ”My F i r s t C program \n” ) ;
}
File → Save
This will open a window with the Save File As option. Enter the desired name of the program, like
first.c in the New File Name and click OK.
Now that you have a program written in C programming language which you have named first.c, let
us see how we can run the program to get the results. Compilation is necessary since the computer
obviously does not understand the english language in which you have written the program. Further,
compilation also takes all the header files (see below) like stdio.h, math.h etc and ”attaches” them to
the program you have written. These header files have predefined functions and routines which make
it easier to use the programming language.
Compilation and Executing or running a program is a two step process.
Page 4 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
will compile the program saved in the text file with the name first.c and create an executable output.x
in the same directory. The extension (.c) is meant to indicate the nature of this file. You can do an ls
to check if output.x is created.
Now we have created an executable file, in our case output.x which can be executed or run by the
machine.
In order to run the executable file produced above, type at the command prompt,
$ ./output.x
Note that ./ is necessary before output.x to make sure that the file output.x in the current directory
is executed.
Important: If you dont give a name to the output file (like we have given above (out-
put.x), then the compiler, uses a default name. The default name is a.out. So for instance,
instead of typing as above, if you compiled your first.c program as
$ gcc first.c
then the executable file created will be a.out and NOT output.x. This can be run by
typing the following at the command prompt:
$ ./a.out
If everything goes fine, then it will print “This is my first C program” on the screen and will return
back to the command prompt.
Page 5 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
main ( )
{
float x , y ;
p r i n t f ( ” Supply t h e v a l u e o f x i n r a d i a n s \n” ) ;
s c a n f ( ”%f ” ,&x ) ;
y = sin (x );
p r i n t f ( ”%f , %f \n” , x , y ) ;
}
After keying in the program above and saving it to a text file sinx.c, compile/execute and run your
programme. Before you compile delete file ”a.out” if it already exists. Then type,
The program with the name sinx.c will compile and create an executable sinx in the same directory.
You can do ls to check if sinx is created. In order to run the executable file,
$ ./sinx
If every thing goes well the message in the seventh line will be printed on the screen and the cursor will
begin flashing in the next line prompting you for input. If you now key in the number 3.14159 (thereby
wanting to compute sin(π)) and press < ENTER >, the number 0.000000 will flash on the screen for
an instant; after that the screen returns to the COMMAND PROMPT.
1. The first line in the above program is a comment and it is not part of the program in the sense
that it is not compiled. You can give a comment anywhere in the program using the format
/* (Comment) */. The comments are very useful in making programs more readable and self-
explanatory.
2. The next two lines beginning with #include append two header files to the programme. # is
called a preprocessor directive. stdio stands for standard input output. The stdio.h file governs
the input via scanf statement and output via the printf statement. There are other means of
input/output about which you will learn later. The math.h header file contains declarations of
in-built mathematical functions like sin, cos, etc.
3. The 4th line indicates the start of the main program. Every program must have the function
main(). The open brace signals the beginning of the main program. Correspondingly the last line
contains a close brace which signals the end of the program.
4. The sixth line contains the declaration that the variables x and y are floating point numbers.
Each variable in a C program has to be declared in accordance with its nature such as floating
point, integer, character, long integer, double float etc. An integer variable, i, is declared as int
Page 6 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
i;. All these declarations must be made at the very beginning of the main program.
It is good programming practice to initialize the variables to some nominal values.
5. The seventh statement prints the message ”Supply the value of x in radians” on the screen and
the \n in inverted commas causes the cursor to go to the beginning of the next line. Any matter
appearing within the inverted commas in a printf statement appears as it is on the
screen except for the special characters such as \n, %f, %d, %6.2f etc. The meaning of these
special characters will be explained as we go along. The messages such as the one in line 7 in the
above program are used as prompts.
6. The eighth statement is the input statement requiring you to supply the value of the variable x.
If you have declared x as a floating point number, then it is best to supply x as a decimal
fraction (For example x = 2.0 instead of x = 2). The &x in this statement corresponds to
the pointer to the memory address of the variable x. (If you do not understand the pointers at
this stage, do not worry, it would become clear later. It must be noted, though, that the pointer
feature is a major advantage of C over other advanced languages like Fortran). The next statement
computes sin(x) by making use of math.h and assigns it to the variable y. As can be expected,
the next statement prints the value of two numbers x and y. In this statement %f describes the
nature of the variable viz., x is a floating point variable denoted by f in %f. For integers we would
use %d. The semi-colon at the end of a line is essential; it denotes the end of an executable
statement. To save space, one can write more than one statement in a line, each separated by a
semi-colon but it is not advisable. Proper indentation of programs is crucial as it can often help
us avoid common syntax errors.
Example 1.3- sinxtab.c (Version 1)
Now let us try to run a variant of the previous program which will tabulate the value of sin(x) at
twenty equally spaced values of x lying between x1 and x2 that you would specify. The listing of this
program is as follows:
Page 7 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
1. As in the case of floating point variables, integer variables must also be explicitly declared
in the beginning. It is possible to assign values to integers in the declaration statement itself,
as is done for n in the above program.
2. The atan function evaluates arc tangent (arctan) of its argument.
3. The for statement sets a sequence i = 0 to n in unit steps. i++ increments i by one and is identical
to i = i + 1. All statements after the braces on the next line until the closing braces (in the 14th
line) are executed for every value of i in the sequence defined above. This procedure is called
looping.
4. %6.2f is a format symbol which stands for a floating point number of total field of six spaces with
two places after decimal.
5. The space between %6.2f and the second %6.2f is to separate the numbers adequately when
printed on the screen.
Change of Type (typecasts)
Very often one has to change a variable from one type to another; for example int to float, float to
double, float to int etc. Thus float(n) OR (float)(n) converts the value of n into float type.
Note that whereas n remains of the type int and its value in the above example remains
20, the value of (float)(n) becomes 20.0. Similarly, the statement (int)(x*y+0.1*exp(x)) will, for
example, evaluate the expression (x*y+0.1*exp(x)), which is of the type float and then convert it into
type int. Thus we may have
i = (int)(x*y+0.1*exp(x));
Of course i must already have been declared as an int variable. If the expression on the RHS evaluates
to 3.7, the integer i will have the value 3; if it evaluates to -3.8, i will have the value -3.
An alternative program
Page 8 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
x2 = x2 ∗ p i ;
dx = ( x2−x1 ) / ( f l o a t ) ( n ) ;
f o r ( x=x1 ; x <= x2 ; x+ = dx )
{ y = sin (x );
p r i n t f ( ” %6.2 f %6.2 f \n” , x , y ) ;
}
}
to set up the sequence: x=x1, x1+dx, ...x2 (=x1+n*dx). Here x+=dx is an abbreviation for x = x
+ dx. Similarly x*=(expression) implies x = x * expression, etc. Note that the integer variable i is no
longer required in the program; therefore it may be deleted from the declaration statement int i, n =
20. Rerun the program. The results of these two programs should be identical.
1.4 FUNCTIONS
You have already used library functions like sin, cos, atan, etc. To be able to appreciate the modular
structure of C language, it is essential for you to learn how to make your own function routine which
performs a specific computation when called by the main program. In this manner, the main program
will simply become a calling program and call various functions to perform their respective computa-
tions.
Example 1.4
/∗
F u n c t i o n w i t h name f o f t y p e f l o a t and has one
argument which i s o f t y p e f l o a t
∗/
float f ( float x)
{
float z ;
z = x∗x∗x + s i n ( x ) ∗ l o g ( x ) ;
return ( z ) ;
}
/∗Main f u n c t i o n which w i l l c a l l a b o v e f u n c t i o n ∗/
Page 9 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
main ( )
{
float y ;
f l o a t f ( f l o a t x ) ; /∗ f u n c t i o n d e c l a r a t i o n ∗/
p r i n t f ( ” Supply t h e v a l u e f o r which you want t o e v a l u a t e t h e f u n c t i o n \n” ) ;
s c a n f ( ”%f ” ,&y ) ;
p r i n t f ( ” %6.2 f %6.2 f \n” , y , f ( y ) ) ;
}
This method of writing a function returns to the calling program the value of the function for the
specified argument x. Note that x is a floating point variable as also is the return type of the function
f. A float function returns a floating point number. Alternatively, we could define the expression as
follows:
float f(float x)
{
return (x*x*x + sin(x)*log(x));
}
Note that in this case we did not use the local variable z to first calculate the result and then return
the value.
(ii) For such one-line functions you can also write preprocessor directives:
There must be
(a) A space between #define and f(x) and between f(x) and the beginning of the function.
(b) The argument of the function should always be within brackets wherever it appears
and the entire function should also be enclosed within brackets.
The return type of inline functions (int/float etc.) is decided by the argument x. The program may
look like the following:
/∗Main f u n c t i o n which w i l l c a l l a b o v e f u n c t i o n ∗/
main ( )
{
float y ;
p r i n t f ( ” Supply t h e v a l u e f o r which you want t o e v a l u a t e t h e f u n c t i o n \n” ) ;
Page 10 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
s c a n f ( ”%f ” , &y ) ;
p r i n t f ( ” %6.2 f %6.2 f \n” , y , f ( y ) ) ;
}
Example 1.5
return ( y ) ;
}
/∗
F u n c t i o n w i t h name f o f t y p e f l o a t and has two
arguments , one o f t y p e f l o a t and o t h e r o f t y p e i n t .
∗/
Page 11 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
if (x < 0.0) {
y = z − 1 +exp ( x ) ;
}
else {
y = z − l o g (1+x ) ;
}
return ( y ) ;
}
/∗ g l o b a l f u n c t i o n d e c l a r a t i o n ∗/
float f ( float x , int n ) ;
/∗ main f u n c t i o n w i l l c a l l f ∗/
main ( )
{
float q ;
int j ;
p r i n t f ( ” Supply t h e v a l u e o f x ( f l o a t ) , n ( i n t ) \ n” ) ;
s c a n f ( ”%f , %d” , &q , &j ) ;
p r i n t f ( ”x=%6.2 f n=%d f =%6.2 f \n” , q , j , f ( q , j ) ) ;
}
Note that here we have defined a function of two variables, one of which, x, is a float variable and the
other, n, is an integer variable. We can, in fact, define in the same way functions of several variables.
The meaning of if ... else statements should be obvious. If there are more than one statements in the
if block, they must be put within . Similarly, for the else block. The general format for if ... else is
the following:
if (condition)
statements ;
else
statements ;
It is possible to have only if block without the else block. You may replace the printf statement by:
p = f(q,j);
If a statement becomes too long, it can be continued in the next line by using \ . A
generalization to functions of many variables is obvious. There are also functions that do not return
any value to the calling program but perform many other tasks. In fact, most of the work in C is done
through such functions. Examples are scanf, printf. Even main () itself is one such function. We
will come across many more examples of such functions in future.
N.B. The function subprogram can be defined either before or after the main program. If it is
defined after the main program, it has to be identified in the main program eg.
Page 12 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
main ( )
{
void f 2 ( f l o a t x , i n t i )
{
float f1 ( float a , float b ) ;
...
...
...
}
main ( )
{
float f1 ( float x , float y ) ;
void f 2 ( f l o a t a , i n t j ) ;
...
...
...
}
Here we have two user-defined functions, f1 of return-type float depending on two input argu-
ments, both of type float, and f2 of return-type void depending on two input arguments, the first of
type float and the second of type void. The prototype of f1 has been declared in f2 as well as in main;
so the function f1 can be used in both of them. However, the prototype of f2 has been declared only in
main and not in f1. Hence, f2 can be used only in main. In general it is best to make a global
prototype declaration:
Page 13 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
immediately after the inclusion of header and other files. Now the two functions can be used in
main as well as in all other programs. Note that the type of the function and the type and the number
of arguments on which the function depends, must match with the prototype declaration; the names of
the arguments can be anything or may not even be included. For example, the function prototype in
the above example could have been
Note on function names: You must always give meaningful function names to user-defined
functions to improve the readability of your code. Try to avoid giving single alphabet
function names.
Wherever the argument in these mathematical functions is a double (double precision variable or
constant), it can be replaced by a float variable or constant but should be avoided. The single pre-
cision counterparts of these functions usually have f appended at the end of their names. For example
the single precision counterpart of sin() function declared in < math.h > is sinf() and so on. The choice
of using single precision (float) or double precision (double) variables for computation depends on the
nature of the problem and accuracy desired.
Page 14 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015
QUESTIONS
PROBLEMS
1. Make a table of the trigonometric functions sin x, cos x, and tan x for 0 ≤ x ≤ π4 .
2. The function f (x, y) of two variables x, y is defined by
Make a table of the function f (x, y) for −1.0 ≤ x, y ≤ 1.0 at intervals of 0.25 for both x and y.
Page 15 of 77;
Chapter 2
2.1 Introduction
GNUPLOT is a command-line program that can generate two- and three-dimensional plots of functions
and data. The program runs on all major computers and operating systems (Linux, UNIX, Windows,
Mac OSX...).
The software is copyrighted but freely distributed (i.e., you don’t have to pay for it). It was originally
intended to function as a software for plotting mathematical functions and data but has outgrown its
credentials. It now supports many non-interactive uses, including web scripting and integration as a
plotting engine for third-party applications like Octave. Gnuplot has been supported and under devel-
opment since 1986.
Gnuplot supports many types of plots in either 2D and 3D. It can draw using lines, points, boxes, con-
tours, vector fields, surfaces, and various associated text. It also supports various specialized plot types.
Gnuplot supports many different types of output: interactive screen terminals (with mouse and
hotkey functionality), direct output to pen plotters or modern printers (including postscript and many
color devices), and output to many types of graphic file formats (eps, fig, jpeg, LaTeX, metafont, pbm,
pdf, png, postscript, svg, ...). Gnuplot is easily configurable to include new devices.
$gnuplot
16
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
On the gnuplot prompt type the following lines one by one (self explanatory)
To turn off the grid, you can “unset grid”, to turn off the xlabel, you can type “set xlabel ’ ’ ”. Type
set at the gnuplot prompt to see all of the options you can turn on and off. To turn on auto-scaling
(without any ‘x’ or ‘y’ labels: default), type “set auto” at gnuplot>.
Note that many of the gnuplot keywords including: using, title, and with can be abbreviated with a
single alphabet as: u, t and w but should be avoided by beginners. Also, each line and point style has
an associated number.
Page 17 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
In order to draw two plots on top of each other, you can replace the last line by
gnuplot>plot [-2*pi : 2*pi] sin(x) t “Sine Wave” with linespoints, cos(x) t “Cosine
Wave” with linespoints
You can exit from gnuplot by typing “exit” (or “quit”) on the gnuplot prompt.
gnuplot>exit
This will create a postscript image file called “plot.eps” of the previous plot. It will be placed in the
same folder in which you are working. You can then use any postscript viewer program like “gv” (or
“evince”) to open your saved graphics file. Remember, the plot will not appear on the screen when you
redirect the terminal type to postscript (first line of the example above), so it may appear as if nothing
has happened. Exit from gnuplot prompt, and then type on the terminal
Page 18 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
Note: once you have finished plotting to a file, you need to set the terminal type back to
x11 (last line of the previous example) so as to view subsequent plots on the screen.
and save it to the gnuplot script file plotscript.p Finally, to execute this script using gnuplot, on
the terminal prompt, type
$gnuplot ./plotscript.p
This will create .eps file plot1.eps. You can directly view this file using gv as in the above example.
Customization
Customization of the axis ranges, axis labels, and plot title, as well as many other features, are
specified using the set command. Specific examples of the set command follow. (The numerical values
used in these examples are arbitrary.) To view your changes type: replot at the gnuplot> prompt at
any time.
Other features which may be customized using the set command are: arrow, border, clip, contour,
Page 19 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
grid, mapping, polar, surface, time, view, and many more. The best way to learn is by reading the
on-line help information, trying the command, and reading the Gnuplot manual.
FUNCTION RETURNS
Bessel, gamma, beta, gamma, and gamma functions are also supported. Many functions can take
complex arguments. Binary and unary operators are also supported. The supported operators in Gnu-
plot are the same as the corresponding operators in the C programming language, except that most
operators accept integer, real, and complex arguments. The ** operator (exponentiation) is supported
as in FORTRAN. Parentheses may be used to change the order of evaluation. The variable names x, y,
and z are used as the default independent variables.
Page 20 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
#n n2 n3
#Comment
1 1 1
2 4 8
3 9 27
4 16 64
5 25 125
6 36 216
7 49 343
8 64 512
9 81 729
10 100 1000
Again start gnuplot in a terminal by writing “gnuplot” on the terminal prompt and type the follow-
ing line
gnuplot> plot “data1.txt” u 1:2 with linespoints,“data1.txt” using 1:3 with linespoints
gnuplot ignores lines starting with # (comment lines.) This will create a plot like this
You can modify these plots again using the x and y-labeling.
You can also save the above plot in a .eps file using the following commands:
Page 21 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
gnuplot>plot “data1.txt” using 1:2 with linespoints, “data1.txt” using 1:3 with lines-
points
gnuplot> set terminal x11
Exit the gnuplot and then type gv plots2.eps on the terminal command prompt to see the output.
/∗ Program f o r e v a l u a t i n g a p e r i o d i c f u n c t i o n & s t o r e v a r i a b l e s i n a
d a t a f i l e , d a t a 2 . t x t ∗/
main ( )
{
float x , y , z ;
FILE ∗ f p=NULL;
f p = f o p e n ( ” data2 . t x t ” , ”w” ) ; /∗ open a f i l e h a n d l e i n w r i t e mode ∗/
f o r ( x= 0 ; x <= 6 ; x+= 0 . 1 )
{
z = x −2.0∗( i n t ) ( x / 2 . 0 ) ;
i f ( z <1.0)
y = z;
else
y = 2.0 − z ;
f p r i n t f ( fp , ”%f \ t %f \n” , x , y ) ; /∗ p r i n t l i n e t o t h e f i l e ∗/
}
f c l o s e ( f p ) ; /∗ c l o s e t h e f i l e h a n d l e ∗/
}
COMMENTS
1. Include the header file <stdio.h> in the main program to use FILE pointers in the main program.
2. Open a new data file in “w” (write) mode with a call to fopen() function. Other modes of opening
Page 22 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
3. In order to store data inside the file, use fprintf() function in the manner we used printf() earlier
for printing the output to the standard output (screen). “\ t” corresponds to a tab space and “\
n” corresponds to new line. Note that fprintf() takes the first argument as the file pointer fp
which points to the file opened earlier in the program.
4. Once data is stored in the file, close the data file using fclose() subroutine.
Now start gnuplot in a terminal and plot the data you saved in the above example:
f (x) = 1 for 0 ≤ x ≤ π
f (x) = −1 for π ≤ x ≤ 2π
f (x + 2π) = f (x)
will generate
Page 23 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
PROBLEMS
x = sin θ, y = A sin(nθ + δ)
for θ in the range 0 ≤ θ ≤ 4π and the following set of parameters
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{ f l o a t x , y1 , y2 , a , n , d ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
You will get a plot like the one shown in Figure 2.6.
Page 24 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
0.8
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-1
-1 -0.5 0 0.5 1
y = x 0≤x<π
= 2π − x π ≤ x < 2π
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t x , y , z , x1 ; i n t n ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
Page 25 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
{y=2∗pi−z ; }
f p r i n t f ( fp , ”%f ,% f \n” , x , y ) ;
}
}
You will get a plot like the one shown in Figure 2.7.
3.5
2.5
1.5
0.5
0
-20 -15 -10 -5 0 5 10 15 20
3. Plot |Θlm (θ)|2 , the square modulus of the orbital wave function for l = 3, m = 0, ±1, ±2, ±3. The
values of |Θlm (θ)| are given by
√
3 14 5 2
Θ3,0 (θ) = cos θ − cos θ
4 3
√
42
sin θ 5 cos2 θ − 1
Θ3,±1 (θ) =
√8
105 2
Θ3,±2 (θ) = sin θ cos2 θ
√ 4
70 3
Θ3,±3 (θ) = sin θ
8
This is a polar plot (r = |Θlm (θ)|2 , θ). You should plot x against y for 0 < θ < 2π where
Page 26 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
/∗ wave f u n c t i o n ∗/
#include <s t d i o . h>
#include <math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t x , x1 , y1 ; i n t n ;
float f ( float x , int n ) ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
p r i n t f ( ” Give v a l u e o f m” ) ;
s c a n f ( ”%d” ,&n ) ;
f o r ( x = 0 . 0 0 0 1 ; x<(2∗ p i ) ; x+=0.005∗ p i )
{
x1=f ( x , n ) ∗ f ( x , n ) ∗ c o s ( x ) ;
y1=f ( x , n ) ∗ f ( x , n ) ∗ s i n ( x ) ;
f p r i n t f ( fp , ”%f ,% f \n” , x1 , y1 ) ;
}
}
f l o a t f ( f l o a t a , i n t m)
{
f l o a t v1 , v2 , v3 ;
i f (m==0)
{ v1 =(3.0∗ s q r t ( 1 4 . 0 ) ) / 4 . 0 ;
v2 = ( 5 . 0 / 3 . 0 ) ∗ pow ( c o s ( a ) ,3) − c o s ( a ) ;
v3=v1 ∗ v2 ;
return ( v3 ) ; }
i f (m==1 | | m==−1)
{ v1=( s q r t ( 4 8 . 0 ) / 8 ) ∗ s i n ( a ) ;
v2 =5.0∗ c o s ( a ) ∗ c o s ( a ) −1;
v3=v1 ∗ v2 ;
return ( v3 ) ; }
i f (m==2 | | m==−2)
{ v1=( s q r t ( 1 0 5 . 0 ) / 4 ) ;
v2=s i n ( a ) ∗ s i n ( a ) ∗ c o s ( a ) ;
v3=v1 ∗ v2 ;
return ( v3 ) ; }
i f (m==3 | | m==−3)
{ v1=( s q r t ( 7 0 . 0 ) / 8 ) ;
v2=s i n ( a ) ∗ s i n ( a ) ∗ s i n ( a ) ;
v3=v1 ∗ v2 ;
return ( v3 ) ; }
}
You will get a plot like the one shown in Figure 2.8.
Page 27 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
Angular Wavefunctions
0.8
m=0
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
-4 -3 -2 -1 0 1 2 3 4
jn (z), n = 0, 1, 2, 3, 4, 5
where
sin(z)
j0 (z) =
z
sin(z) cos(z)
j1 (z) = −
z2 z
and the rest can be obtained by the recurrence relation
jn (z)
jn−1 (z) + jn+1 (z) = (2n + 1)
z
in the range 0 ≤ z ≤ 5 at intervals of 0.01.
#include<s t d i o . h>
#include<math . h>
main ( )
{ f l o a t x , Q0 , Q1 , Q2 ;
int n ;
Page 28 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( x = 0 . 5 0 ; x <=7.0; x+=.01)
{
Q0=s i n ( x ) / x ;
Q1=s i n ( x ) / ( x∗x)− c o s ( x ) / x ;
f p r i n t f ( fp , ”%f \ t %f \ t ” , x , Q0 ) ;
f p r i n t f ( fp , ”%f \ t %f \ t ” , x , Q1 ) ;
f o r ( n=1;n<=5;n++)
{
Q2=(2.0∗ n +1.0)∗Q1/x−Q0 ;
f p r i n t f ( fp , ”%f \ t %f \n” , x , Q2 ) ;
Q0=Q1 ; Q1=Q2 ;
}
}
}
You will get a plot like the one shown in Figure 2.9.
Bessel functions
1
n=0
n=1
n=2
0.8 n=3
n=4
n=5
0.6
0.4
0.2
-0.2
-0.4
0 1 2 3 4 5 6 7
Page 29 of 77;
Chapter 3
3.1 INTRODUCTION
A commonly encountered problem in Physics is the numerical evaluation of mathematical functions like
the exponential, circular, hyperbolic or hypergeometric function for given values of their arguments and
up to some specified accuracy. Most of these functions can be represented by infinite series, infinite
products or continued fractions. Some examples of series representation of a function are given below:
x2 x3
ex = 1 + x + + + ... (3.1)
2! 3!
2 k
∞ − x
x n X 4
Jn (x) = (3.2)
2 k=0
k!(n + k)!
∞
X xj
log(1 + x) = (−1)j−1 (3.3)
j=1
j
We describe below some numerical methods for evaluating both finite and infinite series. [Though
representation of a function in terms of continued fractions and infinite products is also quite useful,
we shall not discuss these here.]
x2 x3 xn
Sn (x) = 1 + x + + + ··· + (3.4)
2! 3! n!
i
Here each term is of the form xi! ; with i = 0, 1, 2, · · · , n. As long as n is a small number, there
is no problem. However, if we wish to find the sum of this series for large n, say n = 20, there is a
serious problem - the computer cannot handle “large” numbers and 20! is a very large” number. We
overcome this problem by not evaluating individual terms of the series. Instead we find the ratio of two
consecutive terms, ti and ti−1 . Suppose this ratio is R. Then ti = Rti−1 . Since R is usually a small
number, it is possible to find all the terms, given the first term t0 , by assigning to i the values 1, 2, 3 · · · .
30
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
By adding these terms we get the required sum. In the specific example of the series Eq.(3.4) above,
i xi−1 ti 2
ti = xi! and ti−1 = (i−1)! . So R = ti−1 = xi . Thus, starting with t0 = 1, we get t1 = x, t2 = x2! , and so
on. The following program can carry out this process:
Program ser1
/∗ program f o r e v a l u a t i n g a f i n i t e s e r i e s ∗/
#include <s t d i o . h>
#include <math . h>
main ( )
{ float x , t , s ;
int n , i ;
p r i n t f ( ” s u p p l y x and t h e number o f terms n \n” ) ;
/∗ i f n=20 , t h e l a s t term i s x / 1 9 ! ∗/
s c a n f ( ”%f ,%d” ,&x,&n ) ;
s =1.0; t =1.0;
/∗ The l a s t l i n e i n i t i a l i z e s sum s and t h e f i r s t term t ∗/
/∗ The f o l l o w i n g l o o p e v a l u a t e s t h e terms and sums them ∗/
f o r ( i = 1 ; i < n ; i ++)
{ t ∗ = x/ i ; s+ = t ; }
p r i n t f ( ” \n” ) ;
p r i n t f ( ”x=%6.2 f , n=%d , sum= %12.5 e ” , x , n , s ) ; }
x
In this program, the statement s+ = t generates the partial sums S2 (x), S3 (x), · · · while t∗ = i
generates the successive terms for i = 1, 2, · · · , n.
Note that the order of the statements t∗ = xi and s+ = t is important. What happens if
they are interchanged? Note also that the initialization s = 1.0 and t = 1.0 must be done
outside the loop over i. What happens if these are done within the loop?
Instead of taking the ratio of ti and ti−1 , we could also take the ratio of ti+1 and ti . In this case,
xi+1 i
ti+1 = (i+1)! and ti = xi! . Thus, t∗ = ti+1
ti
x
= i+1 and i starts from 0 and the first term is t0 = 1. These
two methods are equivalent provided we take care of the initializations of i and t.
printf(”\n”);
shifts the cursor to a new line. This device is used to make the output more readable. If necessary
more than one such statements can be included. In the statement
printf(”x=%6.2f,n=%d,sum= %12.5e”,x,n,s);
the sum s is printed in e-format with a total field of 12 characters with five places after the decimal.
The following examples of numbers expressed in the e-format illustrate its use:
Page 31 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
Number e-format
-.00234 -2.3400e-03
5643.1 5.64310e+03
1.5648 1.56480e+00
All these numbers have a total field of 12 characters, including the signs before the number and after e.
Run the program given above for various values of n (for a fixed x) and for various values of x (for a
fixed n) and study the behaviour of Sn (x) as a function of x and n. Print your results in the form of a
table and also display them in the form of graphs Sn (x) vs. x for fixed n, or Sn (x) vs. n for fixedx.
You would have noticed by now that to obtain the results for a new set of input parameters (x and n
in this case), you have to run the program again. This is an unnecessary irritant and can be avoided by
introducing one more loop outside the loop over i. The outer loop can prompt you to feed new values
of n or x or both. Look at the following program:
Program ser2
/∗ program f o r e v a l u a t i n g a f i n i t e s e r i e s ∗/
#include <s t d i o . h>
#include <math . h>
main ( )
{ float x , t , s ;
int n , i ;
p r i n t f ( ” \n s u p p l y x and t h e number o f terms n \n” ) ;
s c a n f ( ”%f ,%d” ,&x,&n ) ;
do
{ s = 0 . 0 ; t =1;
f o r ( i =2; i<=n ; i ++)
{ s+=t ; t∗=x / ( i −1); }
p r i n t f ( ” \nx=%6.2 f n=%d sum= %12.5 e ” , x , n , s ) ;
p r i n t f ( ” \n\n” ) ;
p r i n t f ( ” e n t e r x and n ( n −ve t o break t h e l o o p ) \ n” ) ;
s c a n f ( ”%f ,%d” ,&x,&n ) ;
}
while ( n >0); }
Notice that after this program has run for one set of values of x and n and printed the result, it
prompts you to input new values of x and n. It also tells you that by feeding zero or a negative value for
n you can break the loop and come out of it, since the loop will run only for positive values of n. This
is indicated by the while statement at the end of the do loop. Remember that in C or C++, a do
loop must end with a while condition. Also note that while the initialization s = 0.0; t = 1.0; is
outside the for loop, it is inside the do loop. WHY?
If you wish to plot Sn (x) against n (for a fixed x) you can do so by modifying the above program. Store
variables of your program in a .txt file and then use gnuplot to plot that data. But now you will have
to change n regularly. A better way to plot would be to replace the do-while loop by a for loop in
which n changes by a fixed step.
Page 32 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
x 2 x3
+ + ···
Sn (x) = 1 + x + (3.5)
2! 3!
This is the infinite series representation of the exponential function. The program for the evaluation
of this series will only be a slight modification of the ser1 program on the following lines:
main ( )
{ float x , t , s , acc ;
int i ;
...
...
do
{ ...
i=i +1;
...
... }
while ( f a b s ( t / s )> a c c ) ;
/∗ acc i s t h e d e s i r e d a c c u r a c y s a y 1 . 0 e−5 ∗/
p r i n t f ( ”x=%6.2 f , sum= %12.5 e , expx =%12.5 e ” , x , s , exp ( x ) ) ; }
As in the example of finite series, modify your program by introducing an outer loop to run the
program for various values of x. Alternatively, make a table of values of x against exp(x) or draw a
graph of x against exp(x) using the series representation.
Page 33 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
The above program is the prototype of all such programs to evaluate convergent infinite series.
3.4 Questions
1. Why do we evaluate ratio of consecutive terms, and not individual terms while evaluating sum of
a series?
2. What is the necessary condition to sum up an infinite series?
3. What are the possible ways of defining accuracy while summing up an infinite series? Which one
is better and why?
4. Which loop do you generally use for summing up a finite series? What about infinite series?
5. Compare the behaviour of some simple functions at large x (like Prob.2 and Prob.3) and explain
the discrepancy.
3.5 Problems
1. Write a program to evaluate the sum up to 20 terms of the series
1 1 1
1+ 2
+ 3 + 4 + ···
x x x
for a given x(0 ≤ x ≤ 2), and compare your result with the analytic sum of the series.
#include<s t d i o . h>
#include<math . h>
#include<s t d l i b . h>
main ( )
{ float x , t , s , f ; int n ;
f o r ( x = 1 . 1 ; x <=5.0; x+=0.25)
{
t =1.0/( x∗x ) ; s=t ; n=1;
do
{ t ∗=1/x ;
s+=t ;
n+=1;
Page 34 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
}
while ( n<=18);
s=1+s ;
f =(x∗x−x +1)/( x ∗ ( x − 1 ) ) ;
/∗ t h e a b o v e f i s t h e a n a l y t i c a l r e s u l t
f o r t h e INFINITE s e r i e s ∗/
p r i n t f ( ”No . o f Terms=%d x=%f sum=%f a n a l y t i c=%f \n” , n+1,x , s , f ) ;
}
}
x2 x4
cos(x) = 1 − + − ···
2! 4!
accurate to four significant places. Plot cos(x) vs x in the range 0 ≤ x ≤ π.
/∗ c o s ( x ) s e r i e s ∗/
/∗ ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ∗/
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t x , t , sum , z , a c c = 0 . 0 0 0 0 1 ; i n t i ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( x =0;x<=(6∗ p i ) ; x=x +0.02)
{
sum = 1 . 0 ; t = 1 . 0 ; i =1;
do
{ t∗=−x∗x / ( 2 . 0 ∗ i ∗ ( 2 . 0 ∗ i − 1 . 0 ) ) ;
sum+=t ;
i ++;
}
while ( f a b s ( t /sum)> a c c ) ;
f p r i n t f ( fp , ”%f ,% f \n” , x , sum ) ;
}
f p r i n t f ( fp , ” \n ”EXACT RESULTS\n” ) ;
f o r ( x = 0 . 0 ; x<=(6∗ p i ) ; x=x +0.02)
{ f p r i n t f ( fp , ”%f ,% f \n” , x , c o s ( x ) ) ; }
}
Page 35 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
Cos(x) series
2
series
exact
1.5
0.5
-0.5
-1
-1.5
0 2 4 6 8 10 12 14 16 18 20
3. Write a program to evaluate Jn (x) to an accuracy of four significant figures using the following
series expansion:
2 k
∞ (−1)k
x n X x
4
Jn (x) =
2 k=0
k!(n + k)!
Plot Jn (x) against x for 0 ≤ x ≤ 10 and n = 0, 1, 2. Compare with the known behaviour of these
functions and explain the discrepancy at large x.
#include<s t d i o . h>
#include<math . h>
int f a c t ( int n)
{
i n t i , f a c =1;
i f ( n==0) f a c =1;
else
{ f o r ( i =1; i<=n ; i ++)
{ f a c ∗= i ; } }
return ( f a c ) ;
}
main ( )
{
Page 36 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
f l o a t x , sum , t , a c c = 0 . 0 0 1 ; i n t k , n ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( n=0;n<=2;n++)
{
f o r ( x =0;x <=15.0; x=x +.1)
{k =0;sum=1.0/ f a c t ( n ) ; t=sum ;
do
{k++;
t∗=−x∗x / ( 4 . 0 ∗ k ∗ ( n+k ) ) ;
sum+=t ;
}
while ( f a b s ( t /sum)> a c c ) ;
sum∗=pow ( ( x / 2 ) , n ) ;
f p r i n t f ( fp , ”%d\ t %f \ t %f \ t %f \n” , n , f a c t ( n ) , x , sum ) ;
}
}
}
0.8
0.6
0.4
0.2
-0.2
-0.4
-0.6
0 2 4 6 8 10 12 14 16
Page 37 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
/∗ chap 2 p r o b 4 ∗/
#include<s t d i o . h>
#include<math . h>
#include<s t d l i b . h>
#define p i 3 . 1 4 1 5 9
main ( )
{ f l o a t z , t , s , a c c =0.0001 , s 1 ; i n t i ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( z = 0 . 0 ; z <=1.0; z +=0.05)
{
t=z ; s=t ; i =1;
do
{ t∗=−p i ∗ p i ∗pow ( z , 4 ) / ( 4 ∗ i +1);
s+=t ;
i +=1;
}
while ( f a b s ( t / ( s+a c c )) > a c c ) ;
s 1=c o s ( p i ∗ z ∗ z / 2 ) ∗ s ;
f p r i n t f ( fp , ”%f \ t %f \n” , z , s 1 ) ;
}
}
z F(z)
0.000000 0.000000
0.050000 0.049999
0.100000 0.099968
0.150000 0.149757
0.200000 0.198976
0.250000 0.246886
0.300000 0.292300
0.350000 0.333530
0.400000 0.368394
0.450000 0.394337
0.500000 0.408678
0.550000 0.409012
0.600000 0.393734
0.650000 0.362603
0.700000 0.317171
0.750000 0.260922
0.800000 0.198957
0.850000 0.137161
0.900000 0.081072
0.950000 0.034735
Page 38 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
∞
X zk
f (z, n) = 1 n−k
k=0,2,4
2n−k k!Γ 2
+ 2
for n = 2 and z in the range 0 ≤ z ≤ 5 . You would require the following relations:
√
1
Γ = π
2
Γ(z + 1) = zΓ(z)
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{ f l o a t z , t , s , a c c = 0 . 0 0 0 1 ; i n t n=2 ,k ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
}
f c l o s e ( fp ) ;
}
Page 39 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
3.5
2.5
1.5
0.5
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
z 3 1.4z 6 1.4.7z 9
f (z) = C 1 + + + + ···
3! 6! 9!
where C = 0.35503, for z in the range −10 ≤ z ≤ 0, at intervals of 0.05.
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{ f l o a t z , t , s , a c c =0.01 , s 1 ; i n t i ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( z = −10.0; z <=2.0; z +=0.05)
{
t = 1 . 0 ; s=t ; i =1;
do
{ t∗=pow ( z , 3 ) / ( 3 ∗ i ∗ ( 3 ∗ i − 1 ) ) ;
s+=t ;
i +=1;
}
while ( f a b s ( t / s )> a c c ) ;
s 1 =0.35503∗ s ;
Page 40 of 77;
FINITE & INFINITE SERIES COMPUTER PROGRAMMING-2015
f p r i n t f ( fp , ”%f \ t %f \n” , z , s 1 ) ;
}
f c l o s e ( fp ) ;
}
0.8
0.6
0.4
0.2
-0.2
-0.4
-10 -8 -6 -4 -2 0 2
Page 41 of 77;
Chapter 4
ROOT FINDING
4.1 INTRODUCTION
Assume that a function f (x) is continuous within the interval (a, b) on the real line. A point x0 within
this interval is said to be a root of the equation f (x) = 0, if f (x0 ) = 0 (see Figure 4.1). We shall discuss
methods of finding all the real roots of a given equation in a specified interval of the variable x.
To find the rough location and the number of roots in the given interval, the most convenient procedure
is to tabulate the function at a sufficiently large number of points within that interval. If the function
has different signs at two points x1 and x2 , then there exists at least one root in the interval (x1 , x2 ).
So, by looking at the changes of sign of the function, one can get an idea of the number and the rough
location of the roots in the required interval.
Once the rough location of a root has been found, our next task is to find the root to a prescribed
accuracy by narrowing down the interval (x1 , x2 ) within which it lies.
Most methods of root finding depend upon what is called the process of iterations. We shall describe
three such methods here. These methods are general and apply to any kind of function. There are
other methods that apply specifically to polynomials and can find even complex roots, but we shall not
discuss them here.
42
ROOT FINDING COMPUTER PROGRAMMING-2015
(xL + xR )
xM =
2
If f (xM ) = 0, then xM is the required root. If f (xM ) is not zero and has the same sign as f (xL )
(Fig.(4.2) then the root must lie between xM and xR . In this case we replace xL by xM and repeat
the process of halving the interval and comparing the signs. In case f (xM ) has a sign opposite
to that of f (xL ) (Fig.(4.3)), then the root must lie between xL and xM . In this case we look for the
root between xL and xM by replacing xR by xM and repeating the entire process. Since the root is
known to lie between xL and xR and this interval (xR xL ) is halved in every successive approximation,
the root can be determined to any desired accuracy. The process of bisecting the interval is repeated
|(xR −xL )|
while |(x R +xL )|
remains greater than the desired accuracy. Since, by definition, at a root the function
must be zero, alternatively the condition for the termination of the loop can depend on the value of
the function itself, i.e., the loop is repeated while |f (xM )| is greater than or equal to some very small
number, say 10−6 .
The bisection method has the advantage that the procedure is guaranteed to converge. The disad-
vantage is that the method is slow. Another disadvantage is that roots lying close to each other may
Page 43 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
be difficult to separate.
As can be seen from the Fig(4.4), x3 provides a better approximation to the actual root x0 . The
process of finding the secant line is repeated with the points x2 and x3 : we use a do while structure
and replace x1 by x2 and x2 by x3 within the loop. The loop is terminated when the magnitude of
the difference between the two successive approximations to the root becomes less than the desired
accuracy. One may, as in the case of the bisection method, put the condition for the termination of the
loop on the value of the function instead of the value of the root.
The secant method usually converges much faster than the bisection method. It has the disadvantage
however that the procedure is not guaranteed to converge.
Page 44 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
Since a is close to a0 , h is expected to be small. If we can find h, then we would know the root exactly,
since a0 = a + h. Now,
h2 00
f (a0 ) = f (a + h) = f (a) + hf 0 (a) + f (a) + · · · (4.2)
2!
If h is sufficiently small, then
h2 00
f (a) hf 0 (a)
2!
We can therefore neglect terms of order higher than h, and since f (a0 ) = 0 (a0 is the actual root),
we get
f (a)
h≈−
f 0 (a)
This provides an estimate of how far the guess value, a, is from the actual root, a0 . In the next
iteration a + h is taken as the new approximation to the root, a0 . This iterative process is continued as
long as | ha | remains greater than the desired accuracy.
It is clear that larger the numerical value of the derivative f 0 (x) in the neighbourhood of the given
root, smaller is the correction h that has to be added to obtain the next approximation to the root.
Newton-Raphson method is therefore particularly convenient when the graph of the function is steep in
the neighbourhood of the given root. On the other hand, if the derivative f 0 (x) is small near the root,
h will be large and computing the root by this method will be difficult (sometimes even impossible).
The main advantage of this method is that only one point is required to start the procedure. Also the
method usually converges very rapidly. However, the disadvantage is that the derivative of the function
has to be calculated, and this may not be easy. It is also not guaranteed that the method will converge
to a root from any arbitrary starting point. In fact some times the method is known to fail even if one
starts reasonably close to the expected root.
4.5 Questions
1. Compare the behaviour of some simple functions at large x (like Prob.2 and Prob.3) and explain
the discrepancy.What do you understand by iterative methods? How do you find out (roughly) if
a given equation has more than one root?
2. How many initial points are needed for the three different methods? How do you select them?
3. What are the advantages and disadvantages of the three methods studied for root-finding?
Page 45 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
4.6 Problems
1. Find the roots, accurate to four significant figures, of the equation
eax − bx2 = 0
in the range −1.0 ≤ x ≤ 1.0 for
i) a= 1.0, b = 5.0
ii)a=-1.5, b=10.0
by three iteration methods, that is Bisection, Secant and Newton-Raphson methods. In each case,
determine the number of iterations necessary to obtain the desired accuracy.
/∗ r o o t s o f exp ( a∗ x)−b ∗ x ∗ x = 0 u s i n g b i s e c ∗/
/∗ u s e s s u b r o u t i m e t o c a l c u l a t e t h e f u n c t i o n ∗/
#include <s t d i o . h>
#include <math . h>
#include<s t d l i b . h>
main ( )
{
f l o a t x l , xm, xr , a , b , d , a c c =0.00001 , z ;
float f ( float x , float a , float b ) ;
Page 46 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
p r i n t f ( ” g i v e v a l u e s o f a , b , x l , xr ” ) ;
s c a n f ( ”%f ,% f ,% f ,% f ” ,&a ,&b,& x l ,& xr ) ;
do
{xm = ( x l+xr ) / 2 . 0 ;
i f ( f (xm)< a c c )
{ p r i n t f ( ”%f ,% f ,% f ,% f ” ,xm, f (xm, a , b ) , a c c ) }
else
{
d=f (xm, a , b ) ∗ f ( x l , a , b ) ;
i f ( d>0)
x l = xm ;
i f ( d< 0 )
xr = xm ;
z = ( x l −xr ) / ( x l+xr ) ; }
while ( f a b s ( z)> a c c ) ;
p r i n t f ( ”%f ,% f ,% f ,% f ” ,xm, f ( x l , a , b ) , f (xm, a , b ) , z ) ;
}
}
/∗ m u l t i p l e r o o t s o f exp ( a∗ x)−b ∗ x ∗ x = 0 u s i n g b i s e c ∗/
#include <s t d i o . h>
#include <math . h>
#include<s t d l i b . h>
main ( )
{
f l o a t x , xmin , xmax , x i n c , t e s t , x l , xm, xr , a , b , d , a c c =0.00001 , z ;
float f ( float x , float a , float b ) ;
p r i n t f ( ” g i v e v a l u e s o f a , b , xmin , xmax , xminc ” ) ;
s c a n f ( ”%f ,% f ,% f ,% f ,% f ” ,&a ,&b,&xmin ,&xmax,& x i n c ) ;
f o r ( x=xmin ; x<=xmax ; x+=x i n c )
{ t e s t=f ( x , a , b ) ∗ f ( x+x i n c , a , b ) ;
i f ( t e s t <0)
{ x l=x ; xr=x+x i n c ;
do
{
xm = ( x l+xr ) / 2 . 0 ;
d=f (xm, a , b ) ∗ f ( x l , a , b ) ;
i f ( d>0)
x l = xm ;
i f ( d< 0 )
xr = xm ;
z = ( x l −xr ) / ( x l+xr ) ;
}
while ( f a b s ( z)> a c c ) ;
Page 47 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
2. Using the series expansion for J0 (x), find its two lowest positive roots to an accuracy of four decimal
places. Note that in this problem you have been asked to find the roots to accuracy of
4 decimal places and not 4 digits.
/∗ F i r s t r o o t o f B e s s e l ( Prob . 2 , Roots ) ∗/
Page 48 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
#include<s t d i o . h>
#include<math . h>
#include<s t d l i b . h>
main ( )
{ f l o a t x , x l , xr , xm, a c c =0.00001 , x i n c =0.5 , z ;
float f ( float x ) ;
f o r ( x = 0 . 0 1 ; x <=4.0; x+=x i n c )
{
i f ( f ( x ) ∗ f ( x+x i n c ) <0)
{ x l=x ; xr=x+x i n c ;
do
{ xm=( x l+xr ) / 2 . 0 ;
/∗
i f ( f a b s ( f (xm))< acc )
{ p r i n t f (”\ n xm=%f f (xm)=% f acc=%f \n ” ,xm , f (xm) , acc ) ; b r e a k ; }
∗/
i f ( f (xm) ∗ f ( x l ) >0)
{ x l=xm; }
i f ( f (xm) ∗ f ( x l ) <0)
{ xr=xm; }
z=f a b s ( ( x l −xr ) / ( x l+xr ) ) ;
p r i n t f ( ”xm=%f f (xm)=%f z=%f a c c=%f \n ” ,xm, f (xm) , z , a c c ) ;
}
while ( z>a c c ) ;
p r i n t f ( ” \ n r o o t=%f f (xm)=%f z=%f a c c=%f \n\n” ,xm, f (xm) , z , a c c ) ;
break ;
}
}
}
float f ( float x)
{ f l o a t s , t , f u n c ; i n t k , n=0;
s = 1 . 0 ; t=s ; n=0;k=0;
do
{k=k+1;
t∗=−pow ( x , 2 ) / ( 4 . 0 ∗ k ∗ ( n+k ) ) ;
s+=t ;
}
while ( f a b s ( t / s ) > 0 . 0 0 0 1 ) ;
return ( s ) ;
}
3. The equation
f (x, y) = 0
defines y as an implicit function of x. . As an example, consider
f (x, y) = x3 + y 3 + xy + 1 = 0
For any given x, this is a cubic equation in y; so y can be found by obtaining the roots (one or
three real roots) of this equation, say by secant method. Plot y as a function of x, for −1.5 ≤
x ≤ 1.5. If for some value of x there are three real roots, (y1 , y2 , y3 ), plot all the three points
(x, y1 ), (x, y2 ), (x, y3 ). You may assume that −2 ≤ y ≤ 2.
Page 49 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
/∗ p l o t s o f f ( x ) = x ∗ x ∗ x + a∗ x + a∗a∗a + 1 ∗/
#include <s t d i o . h>
#include <math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{
float x , a , f ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( a = −1.5; a <=1.5; a=a +0.1)
{
f o r ( x = −2.0; x <=2.0; x=x +0.01)
{ f = x∗x∗x + a ∗x + a ∗ a ∗ a + 1 . 0 ;
f p r i n t f ( fp , ”%f \ t %f \n” , x , f ) ;
}
}}
10
-5
-10
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
/∗ m u l t i p l e r o o t s o f x ∗ x ∗ x + a∗ x + a∗a∗a + 1 = 0 ∗/
#include<s t d i o . h>
Page 50 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
#include<math . h>
main ( )
{
f l o a t a , x , f , g , a i n c =0.1 , x i n c =0.1 , f 1 , f 2 , h , a c c = 0 . 0 0 0 0 1 ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( a = −1.5;a <=1.5; a+=a i n c )
{ f o r ( x = −2.0;x <=2.0; x=x+x i n c )
{ f = x∗x∗x +a ∗x + a ∗ a ∗ a + 1 . 0 ;
g = pow ( ( x+x i n c ) ,3)+ a ∗ ( x+x i n c ) + a ∗ a ∗ a + 1 . 0 ;
i f ( f ∗g <0)
{do
{ f 1 = x∗x∗x + a ∗x + a ∗ a ∗ a + 1 . 0 ;
f 2 = 3 . 0 ∗ x∗x + a ;
h = −f 1 / f 2 ;
x = x+h ;
}
while ( f a b s ( h/x)> a c c ) ;
f p r i n t f ( fp , ”%f \ t %f \ t %f \n” , a , x , f 1 ) ;
}
}
} }
2
roots
1.5
0.5
-0.5
-1
-1.5
-1.5 -1 -0.5 0 0.5 1 1.5
Page 51 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
ψ − sin ψ − ωt = 0
Use the solution to plot the orbit whose radial coordinates are given by
r = a(1 − cos ψ)
cos ψ −
cos θ =
1 − cos ψ
Take ω = 1.0, = 0.8 and a = 2.0. Remember that time t, is only a parameter. The equation has
to be solved for each t in the given interval. For each t, the initial value of ψ can be chosen to be t.
/∗ K e p l e r O r b i t ∗/
/∗ o r b i t f o r omega = 1 . 0 , e= 0 . 8 , and a = 2 ∗/
#include<s t d i o . h>
#include<math . h>
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t f 1 , f 2 , h , x , t , x1 , y1 , r , a c c= 0 . 0 0 0 0 1 , c ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( t = 0 . 0 0 0 1 ; t <=(2.0∗ p i ) ; t=t +0.05)
{x = t ;
do
{ f1 = x − 0.8∗ sin (x) − t ;
f2 = 1.0 − 0.8∗ cos (x ) ;
h = −f 1 / f 2 ;
x = x + h;
}
while ( f a b s ( h/x ) > a c c ) ;
/∗ p r i n t f (”% f ,% f \n ” , t , x ) ; ∗/
}
}
Page 52 of 77;
ROOT FINDING COMPUTER PROGRAMMING-2015
Kepler Orbit
1.5
omega=1,e=0.8 & a=2
0.5
-0.5
-1
-1.5
-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5
Page 53 of 77;
Chapter 5
ORDINARY DIFFERENTIAL
EQUATIONS
5.1 INTRODUCTION
An ordinary differential equation (ODE) is an equation involving derivatives of an unknown function
of one independent variable. ODEs are very important tools in modeling a wide variety of physical
phenomena - oscillating systems, electrical networks, chemical reactions, satellite orbits, etc. Therefore
the solution of ODEs is of great importance in physical sciences, ecology, economics and social sciences.
We examine some of the ways of solving such equations.
dy
y0 = = f (x, y); y(x0 ) = y0
dx
The function f (x, y) and the initial condition (x0 , y0 ) are given. The exact solution of this equation is
a curve in the x − y plane passing through the point (x0 , y0 ) . Solving this differential equation means
finding y(x), the equation representing this curve. In general, we are required to find the solution, y(x),
from x = x0 to some final point, x = xf . For this purpose the region from x0 to xf is divided into n
(x −x )
equal intervals, each of length h = f n 0 . Some approximation procedure is then employed to obtain
yi+1 from (xi , yi ) , where yi = y(xi ),and xi = x0 + ih. Thus, beginning with (x0 , y0 ) one obtains (x1 , y1 )
. The whole process is repeated n times (using for or do while loop) to finally obtain y at x = xn = xf .
Various approximation techniques are available for the numerical solution of differential equations; of
these we shall discuss the Euler’s formula and the Runge-Kutta methods.
54
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
The interval between two consecutive values of x is h, the step size: (h = xi+1 − xi ). Therefore
The Euler method is not recommended for practical use because it is not very accurate when compared
to other, fancier, methods with an equivalent step size. This is because the error term is of order h2 .
Moreover, the method is not very stable. Since it is a very simple method, it is useful to obtain a rough
guess of the expected solution.
Consider however the use of a step like Equation(5.1) to take a trial step to the mid-point of the interval.
Then use the value of both x and y at that mid-point to compute the real step across the whole interval:
k1 = hf (xi , yi )
h k1
k2 = hf (xi + , yi + )
2 2
3
yi+1 = yi + k2 + O(h )
This symmetrization cancels the first order error term, h2 , and makes the method second order, i.e.,
the error term is proportional to h3 . This is called the mid-point or second order Runge-Kutta
method.
The basic philosophy of the Runge-Kutta methods is to find yi+1 by evaluating the derivative of y not
only at the initial point (xi , yi ) , as is done in the Euler method, but also at the end point and some
intermediate points. The increment, ∆y, instead of being just hf (xi , yi ) , is some kind of an average
of all these derivatives. The intermediate points are so chosen that the terms of order h2 , h3 , · · · in the
expression for the error cancel out. So, the error is considerably reduced, being proportional to a much
higher power of h.
The most widely used method in this class of methods is the 4th-order Runge-Kutta method, RK4.
The error in this case is O(h5 ). For the solution of the differential equation
Page 55 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
Page 56 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
5.5 Problems
1. For the differential equation
dy
=x+y ; y(0) = 1,
dx
tabulate y(x) for 1 ≤ x ≤ 5 at intervals of 0.1 for different choices of the step size h (h =
0.01, 0.005, 0.002, 0.0001), along with the analytic solution. Use all the three methods for their
comparative study. Note that though the tabulation is required betweenx = 1 and x = 5
only, the process of solving the equation has to begin from x = 0, since the initial
condition is prescribed at that point. Also notice that tabulation has to be done at
intervals of 0.1 only though the step size, h, is much smaller than that.
/∗ s o l v i n g dy / dx = x+y w i t h
y (0)=1 u s i n g E u l e r ’ s method ∗/
#include<s t d i o . h>
#include<math . h>
main ( )
{
f l o a t x=0,y =1.0 , h =0.05 , s ;
FILE ∗ r e s 1 ;
r e s 1=f o p e n ( ” r e s 1 . dat ” , ”w” ) ;
FILE ∗ r e s 2 ;
r e s 2=f o p e n ( ” r e s 2 . dat ” , ”w” ) ;
do
{
Page 57 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
y = y + h ∗ ( x+y ) ;
x = x + h;
f p r i n t f ( r e s 1 , ”%f %f \n” , x , y ) ;
}
while ( x <=5.0);
f o r ( x=h ; x < 5 . 0 ; x+=h )
{ s =2.0∗ exp ( x)−x − 1 . 0 ;
f p r i n t f ( r e s 2 , ”%f %f \n” , x , s ) ;
}
}
θ00 = − sin θ
The pendulum is released from rest at an angular displacement α i.e. θ(0) = α, θ0 (0) = 0. Use the
RK4 method to solve the equation for α = 0.1, 0.5 and 1.0 and plot θ as a function of time in the
range 0 ≤ t ≤ 8π. Also plot the analytic solution valid in the small θ approximation (sin θ ≈ θ ).
/∗ S o l v i n g y ’ ’ = − s i n y w i t h i n i t i a l c o n d t i o n s y ( x =0) = a l p h a
and y ’ ( x =0) =0 u s i n g RK4 ∗/
#include<s t d i o . h>
#include<math . h>
#define f 1 ( x , u , z ) ( ( z ) )
#define f 2 ( x , u , z ) (( − s i n ( u ) ) ) /∗ t=x , t h e t a = u∗/
main ( )
{
f l o a t h =0.01 , x , u , z , k1 , k2 , k3 , k4 , m1, m2, m3, m4, a ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s 1 . t x t ” , ”w” ) ;
f o r ( a = 0 . 0 ; a <=1.1; a +=0.1)
{
x = 0.0; u = a ; z = 0;
do
{ k1 = h∗ f 1 ( x , u , z ) ;
m1 = h∗ f 2 ( x , u , z ) ;
k2 = h∗ f 1 ( x+h / 2 , u+k1 / 2 , z+m1 / 2 ) ;
m2 = h∗ f 2 ( x+h / 2 , u+k1 / 2 , z+m1 / 2 ) ;
k3 = h∗ f 1 ( x+h / 2 , u+k2 / 2 , z+m2 / 2 ) ;
m3 = h∗ f 2 ( x+h / 2 , u+k2 / 2 , z+m2 / 2 ) ;
k4 = h∗ f 1 ( x+h , u+k3 , z+m3 ) ;
m4 = h∗ f 2 ( x+h , u+k3 , z+m3 ) ;
u = u+(k1 +2.0∗( k2+k3)+k4 ) / 6 . 0 ;
z = z+(m1+2.0∗(m2+m3)+m4 ) / 6 . 0 ;
x = x+h ;
i f (( fabs (a − 0.1) < 0.001) | | ( fabs (a − 0.5) < 0.001) | | ( fabs (a − 1.0) < 0.001))
{ f p r i n t f ( fp , ”%f \ t %f \ t %f \ t %f \n” , x , u , a , a ∗ c o s ( u∗x ) ) ; }
}
while ( x <=1.0);}
}
Page 58 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
Figure 5.1: Pendulum for different initial displacements & harmonic solution
ẋ = γ1 x − γ2 xy, ẏ = −γ3 y + γ4 xy
where x(t) and y(t) represent respectively the prey and predator populations as functions of time.
The term γ1 x tells us that the prey population grows in proportion to its own population while
−γ2 xy says that it decreases as a result of encounters with the predators. The second equation says
that the predator population decreases in proportion to its own population (to model competition
for food between its members) and increases as a result of encounters with the prey (by providing
food for the predators). Solve these equations for γ1 = 0.25, γ2 = γ4 = 0.01 and γ3 = 1 with the
initial values x(0) = 100 and successively y(0) = 5, 10, 15, 20, and 25. Plot x(t) versus y(t) for
0 ≤ t ≤ 20.
/∗ s o l v i n g dx / d t = g1 ∗ x + g2 ∗ x ∗ y
and dy / d t = −g3 ∗ y + g4 ∗ x ∗ y RK4 method ∗/
#include<s t d i o . h>
#include<math . h>
#define f 1 ( t , x , y ) ( g1 ∗ ( x ) − g2 ∗ ( x ) ∗ ( y ) )
#define f 2 ( t , x , y ) (−g3 ∗ ( y ) + g4 ∗ ( x ) ∗ ( y ) )
main ( )
{
f l o a t g1 =0.25 , g2 =0.01 , g3 =1.0 , g4 =0.01 , h = 0 . 0 0 1 ;
f l o a t t , x , y , k1 , k2 , k3 , k4 , m1, m2, m3, m4, y1 ;
Page 59 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
f o r ( y1 = 5 . 0 ; y1 <=25.0; y1 +=5.0)
{ t = 0 . 0 ; x = 1 0 0 . 0 ; y = y1 ;
do
{ k1 = h∗ f 1 ( t , x , y ) ;
m1 = h∗ f 2 ( t , x , y ) ;
k2 = h∗ f 1 ( t+h / 2 , x+k1 / 2 , y+m1 / 2 ) ;
m2 = h∗ f 2 ( t+h / 2 , x+k1 / 2 , y+m1 / 2 ) ;
k3 = h∗ f 1 ( t+h / 2 , x+k2 / 2 , y+m2 / 2 ) ;
m3 = h∗ f 2 ( t+h / 2 , x+k2 / 2 , y+m2 / 2 ) ;
k4 = h∗ f 1 ( t+h , x+k3 , y+m3 ) ;
m4 = h∗ f 2 ( t+h , x+k3 , y+m3 ) ;
x = x+(k1 +2.0∗( k2+k3)+k4 ) / 6 . 0 ;
y = y+(m1+2.0∗(m2+m3)+m4 ) / 6 . 0 ;
t = t+h ;
f p r i n t f ( fp , ”%f \ t %f \n” , x , y ) ;
}
while ( t <=20.0);}
}
70
60
50
40
30
20
10
0
40 60 80 100 120 140 160 180
d2 y dy
f1 (x) 2 + f2 (x) + f3 (x)y = f4 (x)
dx dx
Page 60 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
with
y(0) = 0 at x = 0
dy
= 1 at x = 0
dx
where
10
X (−1)i x2i
f4 (x) =
i=0
(2i + 1)!
main ( )
{ float f ( float x ) ;
f l o a t h =0.01 , x , u , z , k1 , k2 , k3 , k4 , m1, m2, m3, m4 ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
x = 0.0; u = 0.0; z = 1.0;
do
{ k1 = h∗ f 1 ( x , u , z ) ;
m1 = h∗ f 2 ( x , u , z ) ;
k2 = h∗ f 1 ( x+h / 2 , u+k1 / 2 , z+m1 / 2 ) ;
m2 = h∗ f 2 ( x+h / 2 , u+k1 / 2 , z+m1 / 2 ) ;
k3 = h∗ f 1 ( x+h / 2 , u+k2 / 2 , z+m2 / 2 ) ;
m3 = h∗ f 2 ( x+h / 2 , u+k2 / 2 , z+m2 / 2 ) ;
k4 = h∗ f 1 ( x+h , u+k3 , z+m3 ) ;
m4 = h∗ f 2 ( x+h , u+k3 , z+m3 ) ;
u = u+(k1 +2.0∗( k2+k3)+k4 ) / 6 . 0 ;
z = z+(m1+2.0∗(m2+m3)+m4 ) / 6 . 0 ;
x = x+h ;
f p r i n t f ( fp , ”%f \ t %f \n” , x , u ) ;
}
while ( x <=1.0);
}
float f ( float x)
{
f l o a t f =1.0 ,R= 1 . 0 ;
int k ;
f o r ( k=1;k<=10;k+=1)
{ R = − R∗x∗x / ( 2 . 0 ∗ k ∗ ( 2 . 0 ∗ k + 1 . 0 ) ) ;
f = f + R;
Page 61 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
}
return ( f ) ;
}
0.8
Plot of given differential equation
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
5. Do numerical integration on the following differential equations (Lorenz equations) with integration
step size, ∆t = 0.01:
dx dy dz 8
= −10(x − y) = αx − xz − y = xy − z
dt dt dt 3
Plot the trajectories (after removing transients)
a) in x − y, x − z, y − z planes, and
b) in x − t, y − t, z − t planes,
Page 62 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
/∗ s o l v i n g dx / d t = s ( y−z ) ,
dy / d t = x ( r−z )−y ,
and dz / d t = xy−b z u s i n g RK4 method ∗/
#include<s t d i o . h>
#include<math . h>
#define f 1 ( t , x , y , z ) ( s ∗ ( ( y) −(x ) ) )
#define f 2 ( t , x , y , z ) ( ( x ) ∗ ( r −(z )) −( y ) )
#define f 3 ( t , x , y , z ) ( ( x ) ∗ ( y)−b ∗ ( z ) )
main ( )
{
f l o a t s =10.0 , r =5.0 , b = 8 . 0 / 3 . 0 , h = 0 . 0 1 ;
f l o a t t , x , y , z , k1 , k2 , k3 , k4 , m1, m2, m3, m4, n1 , n2 , n3 , n4 , p ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
t = 0.0; x = 1.00; y = 2.0; z = 1.0;
do
{ k1 = h∗ f 1 ( t , x , y , z ) ;
m1 = h∗ f 2 ( t , x , y , z ) ;
n1 = h∗ f 3 ( t , x , y , z ) ;
k2 = h∗ f 1 ( t+h / 2 , x+k1 / 2 , y+m1/ 2 , z+n1 / 2 ) ;
m2 = h∗ f 2 ( t+h / 2 , x+k1 / 2 , y+m1/ 2 , z+n1 / 2 ) ;
n2 = h∗ f 3 ( t+h / 2 , x+k1 / 2 , y+m1/ 2 , z+n1 / 2 ) ;
k3 = h∗ f 1 ( t+h / 2 , x+k2 / 2 , y+m2/ 2 , z+n2 / 2 ) ;
m3 = h∗ f 2 ( t+h / 2 , x+k2 / 2 , y+m2/ 2 , z+n2 / 2 ) ;
n3 = h∗ f 3 ( t+h / 2 , x+k2 / 2 , y+m2/ 2 , z+n2 / 2 ) ;
k4 = h∗ f 1 ( t+h , x+k3 , y+m3, z+n3 ) ;
m4 = h∗ f 2 ( t+h , x+k3 , y+m3, z+n3 ) ;
n4 = h∗ f 3 ( t+h , x+k3 , y+m3, z+n3 ) ;
x = x+(k1 +2.0∗( k2+k3)+k4 ) / 6 . 0 ;
y = y+(m1+2.0∗(m2+m3)+m4 ) / 6 . 0 ;
z = z+(n1 +2.0∗( n2+n3)+n4 ) / 6 . 0 ;
t = t+h ;
p = ( x+y ) / 2 . 0 ;
f p r i n t f ( fp , ”%f \ t %f \ t %f \ t %f \n” , t , x , y , z ) ;
}
while ( t <=20.0);
}
The plots of x vs t for the various values of α will look like the Figures below.
Page 63 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
3.5
2.5
1.5
1
0 2 4 6 8 10 12 14 16 18 20 22
20
10
-10
-20
-30
0 2 4 6 8 10 12 14 16 18 20 22
Page 64 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
40
30
20
10
-10
-20
-30
-40
0 2 4 6 8 10 12 14 16 18 20 22
/∗ s o l v i n g x ( n+1)= a l p h a x n − a l p h a x n ˆ2 ∗/
#include<s t d i o . h>
#include<math . h>
main ( )
{
float a , x ; int i ;
FILE ∗ f p=NULL;
//FILE ∗ f p 1=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
// f p 1=f o p e n (” r e s 1 . t x t ” ,”w ” ) ;
x =0.9;
f o r ( a = 0 . 0 1 ; a <=4.0; a +=0.05)
{ x =0.9;
Page 65 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015
i =0;
do{
x=a ∗x − a ∗x∗x ;
// f p r i n t f ( fp1 ,”% f \ t %f \n ” , x , a ) ;
i ++;
i f ( i >10000)
{
f p r i n t f ( fp , ”%f \ t %f \ t %d\n” , x , a , i ) ;
}
} while ( i <=10100);
}}
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Page 66 of 77;
Chapter 6
INTEGRATION
6.1 INTRODUCTION
Rb
Evaluation of the definite integral f (x)dx is equivalent to finding the area enclosed by the curves rep-
a
resented by the integrand y(x) , the two ordinates at a and b, and the x-axis (Figure (6.1)). Numerical
integration or quadrature, as it is usually called, aims at calculating this area. To do this the required
area is divided into several strips by erecting ordinates in the interval between a and b . The area of
each strip is found by some suitable approximation and these areas are summed to get the total area.
In some procedures the width of the strips is kept equal, while in others it is unequal. It is of course
assumed that the given integral has a finite value.
We can find the area of each trapezium and add. The result is
67
INTEGRATION COMPUTER PROGRAMMING-2015
Zb
(y0 + y1 ) (y1 + y2 ) (yn−1 + yn )
f (x)dx ∼
= h+ h + ··· + h
2 2 2
a
" n−1
#
h X
= y0 + yn + 2 yi
2 1
" n−1
#
h X
= y(a) + y(b) + 2 y(a + ih) (6.1)
2 1
This formula is called the Trapezoidal rule. It can be shown that for trapezoidal rule the leading
term in the error is
(b − a)h2 y 00 (ξ)
Er = − (6.2)
12
where ξ is some point in the interval (a, b) and y 00 denotes the second derivative of y. Clearly, if
the integrand is a linear function of its argument, all derivatives of order two and higher
vanish, and the Trapezoidal rule produces an exact result.
Zb
h
f (x)dx = [y0 + 4(y1 + y3 + · · · + yn−1 ) + 2(y2 + y4 + · · · + y2n−1 ) + yn ]
3
a
n n
2 2
−1
h X X
= y0 + yn + 4 y2i−1 + 2 y2i
3 i=1 i=1
n n
2 2
−1
h X X
= y(a) + y(b) + 4 y(a + (2i − 1)h) + 2 y(a + 2ih) (6.3)
3 i=1 i=1
It is obvious that in this case n must be even. Simpson’s rule generally yields better results than
the Trapezoidal rule because we use a quadratic rather than a linear approximation to each segment of
the curve. The leading term in the error in Simpson’s rule is proportional to the fourth derivative of
y ; therefore it produces exact results in the case of integrands that are general polynomial of order three.
For both methods, the first step is to choose an appropriate interval width h . Equivalently, one can
choose the number of intervals, n . But for certain exceptions, Simpsons rule yields reasonably accurate
result for n ∼ 50100. In the trapezoidal rule, for same kind of accuracy, a much larger number of points
may be needed. Remember that for Simpson’s rule n has to be even and that the various
terms have to be multiplied by proper weights.
Page 68 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
Zb n
X
f (x)dx ≈ wk f (xk ) (6.4)
1 k=1
and choose the 2n unknowns - the points xk and the weights wk - such that we get an exact result
as long as f (x) is an arbitrary polynomial of order (2n − 1), i.e.,
Z+1 n
X
2 2n−1
wk A0 + A1 xk + A2 x2k + · · · + A2n−1 x2n−1
A0 + A − 1x + A2 x + · · · + A2n−1 x dx = k
−1 k=1
A2 A4 A2n−2
= 2 A0 x + + + ··· +
3 5 (2n − 1)
where the last step is got by explicit integration. Since this relation must be satisfied for arbitrary
A1 , A1 , A2 , · · · , A2n−1 , we obtain the conditions
n
X 2
wk x2i−2
k = i = 1, 2, · · · , n
k=1
2i − 1
n
X
wk x2i−1
k =0 i = 1, 2, · · · , n
k=1
These 2n equations can, in principle, be solved for the 2n unknowns, wk and xk . One can show
that the xk s are the zeros of the Legendre polynomial Pn (x). The weights wk are also related to the
Legendre polynomials.
w1 + w2 = 2
2
w1 x21 + w2 x22 =
3
w1 x1 + w2 x2 = w1 x31 + w2 x32 = 0
with the solution
w1 = w2 = 1
1
x1 = −x2 = √
3
Page 69 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
The approximation
Zb n
X
f (x)dx ≈ wk f (xk )
1 k=1
with the wk ’s and xk ’s determined as above, is called the Gauss-Legendre quadrature, or simply
Gauss quadrature. As the above analysis shows, this procedure should, in general, prove to be supe-
rior to the Trapezoidal rule or Simpson’s rule. The accuracy of the quadrature improves with increasing
. For most functions and for reasonable accuracy, n = 8 is usually quite adequate. The weights wk
and the points xk are given in standard tables for different values of n . Since xk are symmetric about
x = 0, and the weight factors for the two symmetric zeroes are the same, only the positive values of xk
and the corresponding wk need be read in. The negative values of xk are generated from the positive
values. Thus, all one has to do in practice is to decide upon the value of n, read wk and xk from the
tables and evaluate the series on the RHS of Equation(6.4) which gives the value of the integral on the
LHS. Actually, one should perform the calculation for two values of n, say n = 6, and n = 8 and verify
that the two results agree with each other to the desired accuracy. If the agreement is not satisfactory,
a still higher value of n should be used.
The Gauss quadrature as described above applies only when the limits of integration are (−1, 1).
However, an integral with any finite limits (a, b) can be transformed into one with limits (−1, 1) by a
linear transformation of the independent variable:
Zb Z+1
(b − a) (a + b) (b − a)
I= f (z)dz = f( + x)dx
2 2 2
a −1
where
(a + b) (b − a)
z= + x
2 2
Now applying Gauss quadrature we have the result
n
(b − a) X
I= wi f (zi )
2 i=1
where
(a + b) (b − a)
zi = + xi (6.5)
2 2
The xi ’s and wi ’s are obtained from the tables, the zi ’s are then obtained from xi ’s with the help
of Equation(6.5), the given function f (z) is evaluated at these points and the series on the RHS is
evaluated. Finally the result is multiplied by the factor (b−a)
2
to obtain the value of the integral.
Page 70 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
is known as Gauss-Laguerre quadrature. [note that only f (xk ) appears on the RHS in
the sum and not the full integrand e−xk f (xk ) .] The weights wk and the points xk are obtained
from considerations similar to the case of Gauss Quadrature, viz., that the series on the RHS give exact
results whenever is an arbitrary polynomial of order 2n − 1 . In the present case it turns out that the
xk ’s are zeros of the Laguerre polynomial Ln (x) and wK ’s are also related to the Laguerre polynomials.
Again, wk and xk are available in standard tables.
where xk ’s are the zeros of the Hermite polynomial Hn (x) and the weights wk ’s are also related to
the Hermite polynomials. As in the case of the Legendre polynomials, zeros of the Hermite polynomials
are also placed symmetrically about the origin and therefore the number of data points to be fed to the
computer is halved. The xk s and wk s can be read from the tables.
To make such a file for various values of n one uses the switch statement. This statement tests a
single expression and provides different actions depending on its value. The general form of switch is
as follows:
switch ( n )
{
case 1 :
[ statements ]
break ;
case 2 :
[ statements ]
break ;
...
}
Page 71 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
In this example, if n has value 1, the statements following case 1 : are executed (notice the colon
after 1). If n = 2, then statements following case 2 : are executed and so on. So for different val-
ues of n, we can store zeros and weights in the form of arrays and can call them for the desired value of n.
The three files, GAUSS.C, LAGUERRE.C and HERMITE.C have been provided to you in the
include directory. These contain the data for the Gauss-Legendre, Gauss-Laguerre and Gauss-Hermite
quadratures respectively, for some specific values of n: n = 4,6,8,10,12 for Gauss-Legendre and Gauss-
Hermite; n = 4,6,8,10 for Gauss-Laguerre quadrature. The required file must be included only
after n has been defined. For example, if you wish to find an integral using the Gauss quadrature
for n = 8 then your program may be
main ( )
{...
n = 8;
#i n c l u d e <GAUSS. C>
Note: For all three quadratures, the points xi and the weights wi are represented by
arrays x[i] and w[i] respectively. Hence, in your program also you are obliged to use the
same variable names, x[i] and w[i]. Also note that the index, i, begins with zero and not 1.
Further, due to the symmetry of the zeros and weights for the Gauss and Gauss-Hermite
quadratures, the index runs from 0 to ( n2 − 1): only half the values are supplied in the
table, the other half have to be generated by symmetry. For Laguerre quadrature the
index runs from 0 to (n − 1).
6.4 Problems
1. Integrate the following functions to an accuracy of 1 in 105 for given limits a and b:
a)
arctan x
, a = 5, b = 10 (0.142208)
x2
Use Trapezoidal, Simpson & Gauss Quadrature.
b)
2
e−x log(1 + x2 ), a = −∞, b = +∞ (0.587080, 10 pt)
c)
e−x log(1 + x), a = 0, b = +∞ (0.596387, 8 pt)
#include<s t d i o . h>
#include<s t d l i b . h>
#include<math . h>
#define f ( x ) ( atan ( x ) / ( ( x ) ∗ ( x ) ) )
main ( )
{ f l o a t c , d , z , z1 , x [ 2 0 ] , w [ 2 0 ] , a =5.0 , b =10.0 , s ;
Page 72 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
i n t i , n=6;
#i n c l u d e <GAUSS. C>
c =0.5∗( b+a ) ; d =0.5∗( b−a ) ;
s =0.0;
f o r ( i =0; i<=n/2 −1; i ++)
{ z=(c+d∗x [ i ] ) ; z1=c−d∗x [ i ] ;
s+=w [ i ] ∗ ( f ( z)+ f ( z1 ) ) ; }
s∗=d ;
p r i n t f ( ” s=%f ” , s ) ;
}
#include<s t d i o . h>
#include<math . h>
#define f ( x ) ( atan ( x ) / ( ( x ) ∗ ( x ) ) )
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t s , a , b , s1 , s2 , h , k ;
i n t i ,N=80000;
a =10.0; b=200.0;
h=(b−a ) /N;
s =0.0;
f o r ( i =1; i<=N−1; i +=2)
{
s+=4∗ f ( a+i ∗h)+2∗ f ( a+( i +1)∗h ) ;
}
s 2=s+f ( a)+ f ( b ) ;
s 2∗=h / 3 ;
p r i n t f ( ” a=%f b=%f i n t e g 2=%f \n” , a , b , s 2 ) ;
}
#include<s t d i o . h>
#include<math . h>
#define f ( x ) ( l o g (1+( x ) ∗ ( x ) ) )
main ( )
{ f l o a t w [ 1 2 ] , x [ 1 2 ] , s1 , s2 , s = 0 . 0 ;
int i , n ;
n=10;
#i n c l u d e <HERMITE. C>
f o r ( i =1; i<=n / 2 ; i ++)
{ s 1=x [ i − 1 ] ;
s 2=−x [ i − 1 ] ;
s +=(( f ( s 1 )+ f ( s 2 ) ) ∗w [ i − 1 ] ) ;
p r i n t f ( ”%f \n” , s ) ;
}
}
#include<s t d i o . h>
Page 73 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
#include<math . h>
#define f ( x ) ( l o g (1+( x ) ) )
main ( )
{ f l o a t w [ 1 4 ] , x [ 1 4 ] , s1 , s = 0 . 0 ;
int i , n ;
n=10;
#include<LAGUERRE. C>
f o r ( i =1; i<=n / 2 ; i ++)
{ s 1=x [ i − 1 ] ;
s+=( f ( s 1 ) ∗w [ i − 1 ] ) ;
p r i n t f ( ”%f \n” , s ) ;
}
}
where Ais the amplitude of oscillations. For small amplitudes it is possible to approximate the
time period to
"
2 #
A
T1 = 2π 1 +
4
Plot T, T1 and the percentage difference between T and T1 as functions of A for 0 < A < π.
#include<s t d i o . h>
#include<math . h>
f l o a t f ( f l o a t A1 , f l o a t z )
{ float p ;
p =1.0/(1.0 − s i n (A1/ 2 ) ∗ s i n (A1/ 2 ) ∗ s i n ( z ) ∗ s i n ( z ) ) ;
return ( p ) ;
}
main ( )
{ f l o a t A1 , a , b , z1 , z2 , s , san , p i =3.14159 , h , P ;
f l o a t f ( f l o a t A1 , f l o a t z ) ;
i n t i , n =1000;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;
Page 74 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
}
}
70
60
50
40
30
20
10
0
0 0.5 1 1.5 2 2.5 3
3. Let R(θ) be the polar coordinates of a particle moving under a central force. Then θ is given as a
function of R by the expression:
ZR
dr
θ(R) = h i 12
2mE
2mV (r) 1
r0 r2 l2
− l2
− r2
Plot the orbit of the particle for V (r) = − kr (inverse square law force). Use Gauss quadrature for
the evaluation of the integral. The upper limit,R is to be varied from r0 to rm , where r0 and rm
Page 75 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
are the two zeroes of the factor in the square brackets. Take m = l = k = 1 and
4. Locate the smallest positive root of the function F (x) , given by:
Zπ
F (x) = cos [xa cos(t)] sin2n+1 tdt
0
Page 76 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015
Z2π
1
Jn (z) = cos(z cos(x))dx
2π
0
The zeroes of the Bessel function in this interval are found to be 2.403629, 3.518313, 8.651533
and 11.788982.
Page 77 of 77;