Teacher Manual

Download as pdf or txt
Download as pdf or txt
You are on page 1of 78

M.Sc.

Computer Programming- 2015


Teacher’s Manual

July 24, 2015


Chapter 1

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.

1.1 GENERAL INTRODUCTION


1.1.1 Operating System
An operating system (OS) is a software that, after being initially loaded into the computer by a boot
program, manages all other programs in a computer and provides programmers/users with an interface
used to access those resources. An OS processes system data and user input, and responds by allocating
and managing tasks and internal system resources as a service to users and programs of the system.
Common contemporary desktop OS are Linux, Mac OS X, Microsoft Windows and Solaris etc.

In this course we shall use the Linux OS.

1.1.1.1 LINUX OPERATING SYSTEM

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.

1.1.3 PROGRAMMING LANGUAGE


The computer language that we shall use is C, distributed by GNU. The reasons for this choice are
threefold:
1. C is a scientific language much like Fortran or C++ or Pascal in structure, but has certain advan-
tages on syntax
2. C has features that allow the programmer to organize programs in a clear, easy, logical way. For
example, C allows meaningful names for variables without any loss of efficiency, yet it gives a
complete freedom of programming style, including flexible ways of making decisions, and a set of
flexible commands for performing tasks repetitively (for, while, do). Also it has a vast library of
mathematical functions that can be harnessed to extensive computations.
3. C is useful not only in scientific computations but is also used by computer scientists who focus
more on machine level programming or web designing. This is obviously a great advantage both
in an academic as well as in a non-academic world.

1.2 A FEW GENERAL LINUX COMMANDS


You need to log on to your accounts using your userid and password. Using mouse, click on (shown
as circled on the figure) Applications → System Tools → Terminal. This will open a terminal as shown
in the figure. The cursor will be on the command prompt [userid@localhost ]$, where userid is your
user name.

Figure 1.1: Screen Shot of opening screen after logging in

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

$cp filename1.c filename2.c


copies one file to another file. Original copy filename1.c will also exist after the operation.

$mv -i filename1.c filename2.c


moves one file to another file. Original copy filename1.c will NOT exist after the operation. The flag -i
ensures that the computer asks (something like Do you want to proceed?) before executing the command.

$rm -i filename1.c
deletes/removes a file. Once again the flag -i ensures a check before the command is executed.

1.3 Introduction to GNU C


On the command prompt of a terminal, type

$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.

1.3.1 CREATE YOUR FIRST C PROGRAMME


This is a three step process:

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

1.3.1.1 EDIT/WRITE YOUR PROGRAMME

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” ) ;
}

After writing this program, click on the taskbar item

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.

1.3.1.2 COMPILE/EXECUTE/RUN YOUR FIRST C PROGRAMME

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.

1. First one has to compile the program to create an executable file


2. Second the compiled executable file has to be run.

Page 4 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015

Let us see how this is done for our program first.c.

In order to compile and create an executable, type

$ g c c <p r o g r a m f i l e n a m e . c> −o <e x e c u t a b l e f i l e n a m e >

So, in our case, for example,

$ gcc first.c –o output.x

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.

1.3.1.3 RUN YOUR FIRST COMPILED C PROGRAMME

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.

CONGRATULATIONS! You are successful in writing, compiling and executing your


first C program.

Page 5 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015

Program 1.2- sinx.c , Version 1


/∗ s i n x . c − This program c a l c u l a t e s s i n ( x ) ∗/
# i n c l u d e <s t d i o . h>
# i n c l u d e <math . h>

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,

$ gcc sinx.c –lm –o sinx

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.

Hints for understanding the above program:

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:

Program 1.3- sinxtab.c , Version 1


/∗ s i n x t a b . c − t a b u l a t e s v a l u e s o f s i n x ∗/
# i n c l u d e <s t d i o . h>
# i n c l u d e <math . h>
main ( )
{
f l o a t x , x1 , x2 , y , dx , p i ;
i n t i , n=20;
p i = 4 . 0 ∗ atan ( 1 . 0 ) ;
p r i n t f ( ” Supply x1 and x2 i n u n i t s o f p i \n” ) ;
s c a n f ( ”%f , %f ” , &x1 , &x2 ) ;
x1 = x1 ∗ p i ;
x2 = x2 ∗ p i ;
dx = ( x2−x1 ) / ( f l o a t ) ( n ) ;
f o r ( i =0; i<= n ; i ++)
{ x = x1 + i ∗ dx ;
y = sin (x );
p r i n t f ( ” %6.2 f %6.2 f \n” , x , y ) ;
}
}

Page 7 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015

Hints for understanding the above program:

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.

Example 1.3 (version 2)

An alternative program

Program 1.3- sinxtab.c , Version 2


/∗ s i n x t a b v e r 2 . c − t a b u l a t e s v a l u e s o f s i n x ∗/
# i n c l u d e <s t d i o . h>
# i n c l u d e <math . h>
main ( )
{
f l o a t x , x1 , x2 , y , dx , p i ;
i n t i , n=20;
p i = 4 . 0 ∗ atan ( 1 . 0 ) ;
p r i n t f ( ” Supply x1 and x2 i n u n i t s o f p i \n” ) ;
s c a n f ( ”%f ,% f ” ,&x1 ,& x2 ) ;
x1 = x1 ∗ p i ;

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

For example if you wish to compute the function

f(x) = x*x*x + sin x * log x

in your program, you can do so in two styles:

(i) Here is one way, type the following:

# include <s t d i o . h>


# include <math . h>

/∗
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:

#define f(x) ((x)*(x)*(x) + sin(x)*log(x))

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:

#include <s t d i o . h>


#include <math . h>
#define f ( x ) ( ( x ) ∗ ( x ) ∗ ( x ) + s i n ( x ) ∗ l o g ( x ) ) /∗ i n l i n e f u n c t i o n ∗/

/∗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

Let us take another example: We wish to evaluate the function

f (x, n) = xn − 1 + ex for x ≤ 0.0


= xn − log(1 + x) for x > 0.0
where n is an integer. Obviously, such a function can not be written in the inline style easily.
However, in the first style it may be written as:
float f ( float x , int n)
{
float y , z =1.0;
int i ;
f o r ( i =1; i<=n ; i ++)
{
z∗=x ;
}
if (x < 0.0)
{
y = z − 1 + exp ( x ) ;
}
else
{
y = z − l o g (1+x ) ;
}

return ( y ) ;
}

and the full program may look like the following:

# include <s t d i o . h>


# include <math . h>

/∗
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 .
∗/

float f ( float x , int n)


{
float y , z =1.0;
int i ;

f o r ( i =1; i<=n ; i ++){


z∗=x ;
}

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.

#include <s t d i o . h>


#include <math . h>

Page 12 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015

main ( )
{

float f ( float x , int n ) ;


.
}

float f ( float x , int n)


{
y=x ;
return ( y ) ;
}

1.5 Function Prototypes


Every user-defined subroutine/function must be declared in the caller function via a function proto-
type. If a function is to be used in many other functions, its prototype must be declared in each one
of them. Alternately, one may make a global declaration of the function prototype by declaring it
outside all functions. For example

float f1 ( float x , float y)


{
...
...
...
}

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:

float f1(float x, float y);


void f2(float x, int n);

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

float f1(float , float );


void f2(float , int );

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.

1.6 Logical Function


&& - logical AND
== - logical Equality
! - logical NOT
!= - not equal to
|| - logical OR

1.7 Mathematical Functions in C- declared in < math.h >


sin (double x); sinh (double x);
cos (double x); cosh (double x);
acos (double x); asin (double x);
tan (double x); tanh (double x);
atan (double x); atan2 (double y, double x);
exp (double x); sqrt (double x);
abs (int x); fabs (double x);
ceil (double x); floor (double x);
labs (double x); fmod (double x, double y);
log (double x); log10 (double x);
pow (double x, double y);

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.

However, if the argument is an integer, it cannot be replaced by a float or double.

Page 14 of 77;
INTRODUCTION COMPUTER PROGRAMMING-2015

QUESTIONS

1. Can we run a program without having function main()?


2. What header file must be included (#include) if you use printf and scanf statements in your
program? What is # in the #include statement?
3. What header file must be used if you use cout, cin and endl statements?
4. Do we need to put a semicolon ; after #include statements?
5. What do you understand by typecasting? How will you convert an int variable into float in an
arithmetic expression?
6. What are the different ways in which a function (other than main() ) can be written and called in
main() program?
7. Where do we use pointer in a simple code?

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

f (x, y) = x2 + y 4 for |x| > |y|


2 2
= x (x + 1) for |x| = |y|
= y 2 + x4 for |x| < |y|

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

Graphics Using GNUPLOT

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.

2.2 PLOTTING USING IN-BUILT FUNCTIONS OF GNUPLOT


We will first demonstrate gnuplot using built-in functions.

2.2.1 Interactive plotting


In the following example, we plot sin(x) between −2π < x < 2π with labels on x and y axis.

Open a terminal and type at the prompt

$gnuplot

You will see a screen like this

16
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015

Figure 2.1: Screen Shot of opening GNUPLOT

On the gnuplot prompt type the following lines one by one (self explanatory)

gnuplot> set xlabel ‘x’


gnuplot> set ylabel ‘sin(x)’
gnuplot> set time
gnuplot> set grid
gnuplot> plot [-2*pi : 2*pi] sin(x) title “Sine Wave” with linespoints

You will get a plot like this

Figure 2.2: Screen Shot of Sin(x)

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

Figure 2.3: Screen Shot of Sin(x) & Cosine(x)

You can exit from gnuplot by typing “exit” (or “quit”) on the gnuplot prompt.

gnuplot>exit

2.3 Saving Plots


To save the above image as an enhanced postscript file with “.eps” extension, instead of displaying it
to the screen, enter the following commands:

gnuplot> set terminal postscript enhanced color solid 22


gnuplot>set output ‘plot.eps’
gnuplot>set xlabel ‘x’
gnuplot>set ylabel ‘y’
gnuplot> plot [-2*pi : 2*pi] sin(x) title “Sine Wave” with linespoints
gnuplot>set term x11

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

$gv plot.eps &

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.

2.4 Plotting using script


You can also plot and save above-mentioned example by writing all the commands in a script. When
this script is executed in gnuplot, each line is executed in sequence beginning from the top. To do that,
open your favourite text editor (gedit, emacs, vi, nedit etc) and then type following lines (exactly the
same which you typed interactively at the gnuplot> prompt):

set terminal postscript enhanced color solid 22


set output ‘plot1.eps’
set xlabel ‘x’
set ylabel ‘y’
plot [-2*pi : 2*pi] sin(x) title “SineWave” with linespoints
set term x11

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.

Create a title: > set title ”Force-Deflection Data”


Put a label on the x-axis: > set xlabel ”Deflection (meters)”
Put a label on the y-axis: > set ylabel ”Force (kN)”
Change the x-axis range: >set xrange [0.001:0.005]
Change the y-axis range: > set yrange [20:500]
Have Gnuplot determine ranges: > set autoscale
Put a label on the plot: > set label ”yield point” at 0.003, 260
Remove all labels: > unset label
Plot using log-axes: > set logscale
Plot using log-axes on y-axis: > unset logscale; set logscale y
Change the tic-marks: > set xtics (0.002,0.004,0.006,0.008)
Return to the default tics: >unset xtics; set xtics auto

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.

The supported functions include:

FUNCTION RETURNS

abs(x) absolute value of x, |x|


acos(x) arc-cosine of x
asin(x) arc-sine of x
atan(x) arc-tangent of x
cos(x) cosine of x, x is in radians.
cosh(x) hyperbolic cosine of x, x is in radians
erf(x) error function of x
exp(x) exponential function of x, base e
inverf(x) inverse error function of x
invnorm(x) inverse normal distribution of x
log(x) log of x, base e
log10(x) log of x, base 10
norm(x) normal Gaussian distribution function
rand(x) pseudo-random number generator
sgn(x) 1 if x > 0, −1 if x < 0, 0 if x = 0
sin(x) sine of x, x is in radians
sinh(x) hyperbolic sine of x, x is in radians
sqrt(x) the square root of x
tan(x) tangent of x, x is in radians
tanh(x) hyperbolic tangent of x, x is in radians

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.

2.5 PLOTTING USING DATA FROM FILE


Create a data file, data1.txt, using your favourite text editor (say, gedit or emacs) and enter the follow-
ing data (there are three columns)

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

Figure 2.4: Screen Shot of Plot using data file

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:

gnuplot>set terminal postscript enhanced color solid 22


gnuplot>set output ‘plots2.eps’

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.

2.6 STORING DATA FILE IN C PROGRAM, TO BE USED WITH


GNUPLOT
In this program, we will first write a C program of a periodic function and then store the data to a
data file data2.txt and plot it using gnuplot. We will learn essential features of storing data file in C
program. Type the following program and give it the name periodic.c

/∗ 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 ∗/

#include <s t d i o . h>


#include <math . h>

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

a file are “a” for append, “r” for read, etc.

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:

gnuplot> plot “data2.txt” using 1:2 with linespoints

2.7 Periodic Function


Very often we need to plot a function which is periodic: y = f (x) = f (x + T ), where T is the period of
the function, over a range which covers many periods. The function may be given in an analytic form
over one period, say, from x = 0 to x = T . For example the periodic function

f (x) = 1 for 0 ≤ x ≤ π
f (x) = −1 for π ≤ x ≤ 2π
f (x + 2π) = f (x)
will generate

Figure 2.5: Screen Shot of Plot of periodic function

To find the value of the function at any arbitrary


 point x, we have to find the “corresponding” point
x
in the range (0, 2π). Now x contains n = int 2π multiples of 2π . Thus x1 = x − 2πn, lies between
(0, 2π), and therefore f (x1) is known. But f (x) = f (x − 2πn) = f (x1) and therefore f (x) is also known.
In the present example if x ≤ π , f (x) = f (x1) = 1; if x > π, f (x) = f (x1) = −1.

Page 23 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015

2.7.1 The functions ceil & floor


AS we have seen above, (int) simply truncates the value of the expression (3.7 to 3 or -3.8 to -3). There
are two functions ceil and floor which have a slightly different purpose. The value of the function
ceil(x) is an integer greater than or equal to x. Similarly, the value of the function floor(x) is an
integer smaller than or equal to x. For example
ceil(2.2) = 3;
floor(2.7) = 2;
ceil(-2.7) = -2;
floor(-2.2) = -3

PROBLEMS

1. Plot the following functions:

x = sin θ, y = A sin(nθ + δ)
for θ in the range 0 ≤ θ ≤ 4π and the following set of parameters

(a) δ = π4 , A = 1 and n = 2, 2.5, 3 (n = 1 gives an ellipse. Check.)


(b) δ = π4 , n = 2 and A = 0.5, 1, 2
(c) n = 2, A = 1 and δ = π4 , π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 , y1 , y2 , a , n , d ;
FILE ∗ f p=NULL;
f p=f o p e n ( ” r e s . t x t ” , ”w” ) ;

p r i n t f ( ”GIVEN VALUES OF A, N AND DELTA” ) ;


s c a n f ( ”%f %f %f ” ,&a ,&n,&d ) ;
f o r ( x =0;x<=4∗p i ; x+=p i / 1 0 0 . 0 )
{ y1 = s i n ( x ) ; y2=a ∗ s i n ( n∗x+d ) ;
f p r i n t f ( fp , ”%f ,% f \n” , y1 , y2 ) ;
}
}

You will get a plot like the one shown in Figure 2.6.

Page 24 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015

Lissajous Figures with A=1,n=2,delta=pi/4


1

0.8

0.6

0.4

0.2

-0.2

-0.4

-0.6

-0.8

-1
-1 -0.5 0 0.5 1

Figure 2.6: Lissajous Figures

2. The periodic function y(x) with period 2π is defined as:

y = x 0≤x<π
= 2π − x π ≤ x < 2π

Plot this function from x = −6π to x = 6π .

#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” ) ;

f o r ( x=−(6∗ p i ) ; x <=(6.0∗ p i ) ; x=x +0.001∗ p i )


{n = ( i n t ) ( x / ( 2 . 0 ∗ p i ) ) ;
x1 = x − 2 . 0 ∗ p i ∗n ;
z=f a b s ( x1 ) ;
i f ( z >=0 && z <p i )
{y=z ; }
i f ( z>= p i && z < 2∗ p i )

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

Figure 2.7: Saw tooth wave

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

x = |Θlm (θ)|2 cos θ

Page 26 of 77;
Graphics Using GNUPLOT COMPUTER PROGRAMMING-2015

y = |Θlm (θ)|2 sin θ

/∗ 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

Figure 2.8: Square Modulus of Orbital Wavefunction

4. Plot on the same graph the functions

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

Figure 2.9: Bessel Functions

Page 29 of 77;
Chapter 3

FINITE & INFINITE SERIES

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.]

3.2 Finite Series


Before looking at the evaluation of infinite series, let us see how to evaluate finite series. We take the
example of the following series with a total of (n + 1) terms:

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.

In the above program the statement:

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

3.3 Infinite Series


Whereas a finite series can always be summed in principle, the sum of an infinite series has a meaning
only if the series is convergent. So it must be ensured that the series under consideration is indeed
convergent before one embarks on its evaluation. For finite series, the number of terms to be summed
is given in advance and a for loop can be used to evaluate the sum. In the case of an infinite series,
obviously an infinite number of terms cannot be summed; rather the series is summed till a desired
accuracy is achieved. One measure of accuracy can be in terms of the number of decimal places up
to which the result is required to be correct. A better method of defining accuracy, and the one we
normally use, is in terms of the number of significant figures up to which the result is required to
be accurate. (Why is this method better?) In either case, the number of terms needed to obtain the
desired accuracy is determined during the process of evaluation itself and is not known in advance. The
simplest procedure, therefore, is to replace the for loop by a do-while loop. The process of
addition of terms continues as long as the magnitude of relative contribution of the term to be added,
i.e., —term/sum— remains larger than the desired accuracy. The strategy for finding successive terms
remains the same.

Let us take the example of the following infinite series:

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:

Program inf ser1


/∗ program f o r e v a l u a t i n g an i n f i n i t e s e r i e s ∗/
#i n c l u d e <s t d i o . h>
#i n c l u d e <math . h>

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 ) ;
}
}

2. Evaluate cos(x) using the series

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 ) ) ; }
}

The plot be like in Figure 3.1

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

Figure 3.1: Cos(x) series sum & exact

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 ) ;
}
}
}

The plot be like in Figure 3.2

Bessel functions series sum


1

0.8

0.6

0.4

0.2

-0.2

-0.4

-0.6
0 2 4 6 8 10 12 14 16

Figure 3.2: Bessel Functions series sum

4. . Evaluate F (z) given by



πz 2 (−1)n π 2n z 4n+1
 X
F (z) = cos
2 n=0
1.5.9. · · · .(4n + 1)
correct to four significant figures, for 0 ≤ z ≤ 1 , at intervals of 0.1.

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 ) ;
}
}

The result is given in Table 3.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

Table 3.1: F (z) vs z

5. Write a program to plot the sum of the following series:

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 o r ( z =0; z <=5.0; z +=0.02)


{
t = 1 . 0 / ( 2 . 0 ∗ s q r t ( p i ) ) ; s=t ; k=2;
do
{ t ∗=( z ∗ z ∗ 2 ∗ ( n−k + 1 . 0 ) ) / ( k ∗ ( k − 1 . 0 ) ) ;
s+=t ;
k+=2;
}
while ( f a b s ( t / s )> a c c ) ;
f p r i n t f ( fp , ”%f \ t %f \n” , z , s ) ;

}
f c l o s e ( fp ) ;
}

The plot would look like Figure 3.3.

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

Figure 3.3: Plot of f (z, n) vs z for n = 2

6. Write a program to plot the following function:

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 ) ;
}

The plot will look as in Figure 3.4

0.8

0.6

0.4

0.2

-0.2

-0.4
-10 -8 -6 -4 -2 0 2

Figure 3.4: Plot of f (z) vs z

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.

Figure 4.1: Root of a function

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

4.2 BISECTION METHOD


Suppose the function f (x) varies as shown in the Figures(4.2 & 4.3) below. Let two points xL and xR
be on opposite sides of the root x0 . In the bisection method, one determines the sign of the function at
the mid-point, xM , of xL and xR :

(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 .

Figure 4.2: Bisection Method

Figure 4.3: Bisection Method

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.

4.3 Secant Method


The Secant method uses a secant line to approximate f (x) in the neighbourhood of a root. In this
method also we need two points to begin with. It differs from the bisection method in that the two
points need not be on the opposite sides of the expected root as long as they are sufficiently close to
it. However it is always better to choose the points to lie on the opposite sides of the root. Let x1
and x2 be the two points close to the root x0 . Consider the two points on the curve y = f (x) having
coordinates (x1 , y1 ) and (x2 , y2 ), where y1 = f (x1 ) and y2 = f (x2 ). The equation of the secant line
through these points is given by
y2 − y1
y − y2 = (x − x2 ) (4.1)
x2 − x1
which intersects the x axis at the point x3 where
x1 y2 − x2 y1
x3 =
y2 − y1

Figure 4.4: Secant Method

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.

4.4 Newton-Raphson Method


This method is based on the Taylor expansion of a function. In geometrical terms it is equivalent to
using the tangent to a curve at a given point rather than the secant between two points to find the root.
Let x = a0 be a root of the equation f (x) = 0, let a be a close approximation to it, and let h = a0 a.

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 ( ax)−b ∗ x ∗ x =0 u s i n g Newton−Raphson ∗/


#include<s t d i o . h>
#include<math . h>
#include<s t d l i b . h>
main ( )
{
float a , b , x , f1 , f2 , h ;
float acc =0.00001; int l s ;
p r i n t f ( ” g i v e v a l u e s o f a , b , x\n” ) ;
s c a n f ( ”%f ,% f ,% f ” ,&a ,&b,&x ) ;
l s =0;
do
{ l s=l s +1;
f 1 = exp ( a ∗x ) − b∗x∗x ;
f 2 = a ∗ exp ( a ∗x ) − 2∗b∗x ;
h = −f 1 / f 2 ;
p r i n t f ( ” s t e p=%d , x=%f , f u n c t i o n=%f , h=%f \n” , l s , x , f 1 , h ) ;
x+=h ; }
while ( f a b s ( h/x ) > a c c ) ;

p r i n t f ( ” a=%f , b=%f , Root=%f \n” , a , b , x ) ;


}

/∗ 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 ) ;
}
}

float f ( float x , float p , float q)


{
float v ;
v=exp ( p∗x)−q∗x∗x ;
return ( v ) ;
}

/∗ 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

p r i n t f ( ”%f ,% f ,% f ,% f \n” ,xm, f ( x l , a , b ) , f (xm, a , b ) , z ) ;


}
}
}
float f ( float x , float p , float q)
{
float v ;
v=exp ( p∗x)−q∗x∗x ;
return ( v ) ;
}

/∗ r o o t s o f exp ( a∗ x)−b ∗ x ∗ x =0 u s i n g s e c a n t method ∗/


#include<s t d i o . h>
#include<math . h>
#include<s t d l i b . h>
#define f ( x , a , b ) ( exp ( ( a ) ∗ ( x )) −( b ) ∗ ( x ) ∗ ( x ) )
main ( )
{
f l o a t a , b , x1 , x2 , x3 , a c c =0.001 , z ; i n t l s =1;
p r i n t f ( ”ENTER a , b , x1 and x2 \n” ) ;
s c a n f ( ”%f ,% f ,% f ,% f ” ,&a ,&b,&x1 ,& x2 ) ;
p r i n t f ( ”%f ,% f ,% f ,% f \n” , a , b , x1 , x2 ) ;
do
{
x3 = ( x1 ∗ f ( x2 , a , b)−x2 ∗ f ( x1 , a , b ) ) / ( f ( x2 , a , b)− f ( x1 , a , b ) ) ;
x1 = x2 ; x2 = x3 ;
p r i n t f ( ” l s=%d x1=%5.2 f x2=%5.2 f x3=%5.2 f f ( x3 )=%6.3 f \n” \
, l s , x1 , x2 , x3 , f ( x3 , a , b ) ) ;
z = f a b s ( x1−x2 ) ;
l s ++;
}
while ( z>a c c ) ;
p r i n t f ( ” \n\n%f %f \n” , x3 , f ( x3 , a , b ) ) ;
}

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 ) ;
}
}}

The plot of the function will look as in Figure 4.5

Plot of the function


15
cubic function of x& y

10

-5

-10
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

Figure 4.5: Plot of the function

/∗ 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 ) ;
}
}
} }

The plot will look as in Figure 4.6

2
roots

1.5

0.5

-0.5

-1

-1.5
-1.5 -1 -0.5 0 0.5 1 1.5

Figure 4.6: Plot of the roots

4. Choosing equally spaced values of t in (0, 2π


ω
), solve the Kepler equation for ψ

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 ) ; ∗/

r = 2.0 ∗ (1.0 − 0.8 ∗ cos (x ) ) ;


c = ( cos (x ) −0.8)/(1.0 − 0.8 ∗ cos (x ) ) ;
x1 = r ∗ c ;
y1 = r ∗ s q r t ( 1 . 0 − c ∗ c ) ;
f p r i n t f ( fp , ”%f \ t %f \n” , x1 , y1 ) ;
f p r i n t f ( fp , ”%f \ t %f \n” , x1 ,−y1 ) ;

}
}

The plot for ω = 1,  = 0.8 and a = 2 will look as in Figure 4.7

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

Figure 4.7: Kepler Orbit

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.

Consider the general first-order differential equation:

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.

5.2 Euler’s Formula


Consider the equation y 0 = f (x, y) . If yi and yi+1 are the solutions at points xi and xi+1 , then we may
obtain yi+1 from the approximation
yi+1 − yi
≈ y 0 ≈ f (xi , yi )
xi+1 − xi

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

yi+1 = yi + hf (xi , yi ) + O(h2 )


or
∆yi = yi+1 − yi ≈ hf (xi , yi ) (5.1)
neglecting O(h2 ), the error term of order h2 . Beginning with (x0 , y0 ) , this gives y1 = (y0 + ∆y0 ) at
the point x1 = (x0 + h). Similarly y2 can be calculated starting with (x1 , y1 ) , and so on. Notice that
this last equation is essentially a Taylor’s expansion of y(x) with only the first two terms retained.

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.

5.3 Runge-Kutta methods


A straightforward extension of this method for better accuracy would be to retain higher order terms
in Taylor’s expansion. However, the calculation of higher order derivatives that this straightforward
extension requires is often cumbersome; as a result, this straightforward extension is rarely employed
in practice.

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

y 0 = f (x, y); y(x0 ) = y0


the RK4 method leads to the following formula:

Page 55 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015

yi+1 ≡ y(xi+1 ) = y(xi + h) = yi + ∆yi


where the increment ∆yi is given by
1
∆yi = (k1 + 2k2 + 2k3 + k4 )
6
where
k1 = hf (x, y)
h k1
k2 = hf (x + , y + )
2 2
h k2
k3 = hf (x + , y + )
2 2
k4 = hf (x + h, y + k3 )
Just as in the Euler’s method, here also beginning with (x0 , y0 ), we obtain y1 = y0 + ∆y, i.e.,
y(x1 ), x1 = x0 + h; then from y1 we obtain y2 and so on. The method is very simple to implement
as a function which calculates yi+1 given (xi , yi ) and h. The accuracy of the result can be checked by
repeating the calculation with a larger number of steps, say 2n, and observing the resulting convergence.

5.4 Simultaneous Equations of First-Order


The Runge-Kutta method (or the Eulers method for that matter) can be extended to a system of any
number of coupled first-order or higher order differential equations. To solve a pair of simultaneous
differential equations:
dy dz
y0 = = f1 (x, y, z); z0 = = f2 (x, y, z)
dx dx
y(x0 ) = y0 z(x0 ) = z0
we evaluate the quantities

k1 = hf1 (xi , yi , zi ) m1 = hf2 (x1 , y1 , zi )


   
h k1 m1 h k1 m1
k2 = hf1 xi + , yi + , zi + m2 = hf2 xi + , yi + , zi +
2 2 2 2 2 2
   
h k2 m2 h k2 m2
k3 = hf1 xi + , yi + , zi + m3 = hf2 xi + , yi + , zi +
2 2 2 2 2 2
k4 = hf1 (xi + h, yi + k3 , zi + m3 ) m4 = hf2 (xi + h, yi + k3 , zi + m3 )
This leads to (at xi+1 = xi + ∆h)
1
yi+1 = yi + (k1 + 2k2 + 2k3 + k4 )
6
1
zi+1 = zi +
(m1 + 2m2 + 2m3 + m4 )
6
Notice that any second-order differential system

Page 56 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015

y 00 = f (x, y, y 0 ) y(x0 ) = α, y 0 (x0 ) = β


can be written as two coupled first order equations by putting y 0 = z

y 0 = z, z 0 = f (x, y, z), y(x0 ) = α, z(x0 ) = β


This is a special case of the system considered above. In general, an nth order ODE can
be converted into a system of n first order ODEs in a like manner.

In case of three simultaneous first order ODE’s:

y 0 = f1 (x, y, z, u), z 0 = f2 (x, y, z, u), u0 = f3 (x, y, z, u)


we introduce one more set of quantities in addition to ki s and mi s, and calculate these in a similar
manner. Extension of the method to any number of coupled ODE’s is obvious.

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 ) ;
}
}

2. The ODE describing the motion of a pendulum is

θ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

The plot will look like as in Figure 5.1

Pendulum for different initial displacements


1
RK4
harmonic solution
0.9

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

3. A simple ”prey-predator” system is modeled by the set of equations

ẋ = γ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);}
}

The plot will look like as in Figure 5.2

Prey-predator system for different initial conditions


80
prey-predator

70

60

50

40

30

20

10

0
40 60 80 100 120 140 160 180

Figure 5.2: Prey-predator system for different initial conditions

4. Solve the following differential equation:

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

f1 (x) = f2 (x) = 1, f3 (x) = 4x


and

10
X (−1)i x2i
f4 (x) =
i=0
(2i + 1)!

and plot the result from x = 0 to x = 1.

/∗ S o l v i n g y ’ ’+ y ’ + 4∗ x ∗ y = Sum( i =1;10) [( −1)ˆ{ i } x ˆ{2 i }/(2 i + 1 ) ! ]


w i t h i n i t i a l c o n d t i o n s y ( x =0) = 0 and y ’ ( x =0) =1 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 ) ( f ( ( x ) ) − 4 . 0 ∗ ( x ) ∗ ( u) −( z ) ) /∗ y =u∗/

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 ) ;
}

The plot will look like as in Figure 5.3

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

Figure 5.3: Plot of given Differential Equation

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,

for the following values of the parameter α:

i)α = 5.0 (fixed point solution)


ii) α = 50.0, 125.0, 200.0 (chaotic motion)
iii) α = 100.0, 150.0, 250.0 (periodic motion)

Page 62 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015

Choose any reasonable initial conditions.

/∗ 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

Lorenz Equations: Fixed Point Solution


4.5
alpha = 5.0

3.5

2.5

1.5

1
0 2 4 6 8 10 12 14 16 18 20 22

Figure 5.4: Lorenz Equations: Fixed Point Solution

Lorenz Equations: Chaotic Motion


30
alpha = 50.0

20

10

-10

-20

-30
0 2 4 6 8 10 12 14 16 18 20 22

Figure 5.5: Lorenz Equations: Chaotic Solution

Page 64 of 77;
ORDINARY DIFFERENTIAL EQUATIONS COMPUTER PROGRAMMING-2015

Lorenz Equations: Periodic Motion


50
alpha = 100.0

40

30

20

10

-10

-20

-30

-40
0 2 4 6 8 10 12 14 16 18 20 22

Figure 5.6: Lorenz Equations: Periodic Solution

6. To plot the bifurcation diagram for the logistic map:


A difference equation is a particular form of recurrence relation which is derived from a differential
equation. Consider a difference equation

xn+1 = αxn (1 − xn ) , x0 ∈ [0, 1]


Here α is a parameter α ∈ [0, 4]. Choose a single initial value x0 of x in the range given for x. This
can be any value EXCLUDING 0, 1 and 0.5. For this value of x0 , solve the difference equation for
various values of α in the range given, taking ∆α = 0.05. Thus you will have 100 values of α. For
the solution to the equation for each α, remove the first 104 iterations since these are transients.
Keep the next 100 iterations for each α and plot a graph between x and α.

/∗ 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);
}}

The plot will look like as in Figure 5.7

Bifurcation Diagram for Logistic Map


1

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

Figure 5.7: Bifurcation Diagram for the Logistic Map

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.

6.2 Methods Based on Intervals of Equal Width


6.2.1 Trapezoidal Rule
Let the range of integration, a to b , be divided into n equal intervals, each of width h = (b−a) n
.
Let the ordinates at the (n + 1) points, x0 = a, x1 = a + h, · · · , xn = an + h = b, be denoted by
y0 = y(x0 ), y1 = y(x1 ), · · · , yn = y(xn ). If we now approximate each segment of the curve between two
successive ordinates by a straight line, each strip becomes a trapezium.

Figure 6.1: Trapezoidal Rule

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.

6.2.2 Simpson’s Rule


A method that gives reasonably accurate results and is quite often used for numerical integration is
Simpson’s rule. If in the previous example, we approximate the function over two adjoining segments
by a quadratic rather than a linear function, we obtain what is known as Simpson’s rule:

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

6.3 METHODS BASED ON INTERVALS OF UNEQUAL WIDTH


6.3.1 GAUSS QUADRATURE
We have noted above that the trapezoidal rule gives exact results for polynomials of order one, whereas
Simpson’s rule gives exact results for polynomials of order three. We can do better than this if we
do not insist on the integrand being evaluated at equal intervals. Suppose we approximate any finite
integral by:

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.

As a simple illustration, for n = 2 one gets the equations

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

We see that x1 , x2 are indeed zeroes of P2 (x).

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

6.3.2 GAUSS-LAGUERRE QUADRATURE


Approximating integrals over the semi-infinite interval (0, ∞), by a sum:
Z∞ n
X
e−x f (x)dx ≈ wk f (xk )
0 k=1

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.

6.3.3 GAUSS-HERMITE QUADRATURE


This quadrature is used to evaluate integrals over the interval (−∞, ∞) using the approximation:
Z∞ n
−x2
X
e y(x)dx ≈ wk y(xk )
−∞ k=1

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.

6.3.4 General Considerations & Programming Tips


Keying in the data every time we wish to evaluate an integral is not a good programming practice. One
should key the data in a file and include this file in the program using the #INCLUDE metastate-
ment. Since Gauss-Legendre, Gauss-Laguerre and Gauss-Hermite quadratures use different sets of zeros
and weights, it is best to make three include files and use them according to the integral to be evaluated.

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 ) ;
}
}

2. The time period of a pendulum is given by the integral


π
Z2
1
T =4 2
dz
A
sin2 x

1 − sin 2
0

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” ) ;

f o r (A1= 0 . 1 ; A1<p i ; A1+=0.1)


{ b=p i / 2 . 0 ; a =0; s =0;h=(b−a ) / n ;

f o r ( i =0; i<=n−1; i +=2)

Page 74 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015

{ s+=4∗ f (A1 , a+i ∗h)+2∗ f (A1 , a+( i +1)∗h ) ;


}

s=s+f (A1 , a)+ f (A1 , b ) ;


s ∗=4∗h / 3 . 0 ;

san =2.0∗ p i ∗(1+(A1 / 4 . 0 ) ∗ ( A1 / 4 . 0 ) ) ;


p r i n t f ( ”%f ,% f ,% f \n” , s , san , A1 ) ;
P=(s−san ) / s ;
f p r i n t f ( fp , ”%f \ t%f \ t %f \ t %f \n” ,A1 , s , san , P ) ;

}
}

The plot will look like as in Figure 6.2

Pendulum Time period vs Amplitude


90
Integrated
Approximation
80 percentage difference

70

60

50

40

30

20

10

0
0 0.5 1 1.5 2 2.5 3

Figure 6.2: Time period of pendulum vs Amplitude

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

i) E = −0.25 ( This gives r0 = 0.6, rm = 3.4 approximately)


ii)E = 0(r0 ≈ 0.5, rm ≈ 5)

4. Locate the smallest positive root of the function F (x) , given by:


F (x) = cos [xa cos(t)] sin2n+1 tdt
0

to an accuracy of 4 significant figures, for n = 1 and a = 1.5.

#include <s t d i o . h>


#include <math . h>
#define g ( x , z ) ( ( c o s ( pow ( x , 1 . 5 ) ∗ c o s ( ( z ) ) ) ) ∗ pow ( s i n ( ( z ) ) , 3 ) )
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t x , xmin=0 ,xmax =10.0 , x i n c =0.5 , t e s t , x l , x1 , x2 , x3 , a c c =0.00001 , z ;
i n t f a c t ( i n t n ) ; f l o a t f ( f l o a t x ) ; i n t l s =0;
f o r ( x=xmin ; x<=xmax ; x+=x i n c )
{ t e s t=f ( x ) ∗ f ( x+x i n c ) ;
i f ( t e s t <0)
{ l s =1; x1=x ; x2=x+x i n c ;
do
{ x3 = ( x1 ∗ f ( x2)−x2 ∗ f ( x1 ) ) / ( f ( x2)− f ( x1 ) ) ;
x1 = x2 ; x2 = x3 ; z = f a b s ( x1−x2 ) ; l s ++;
}
while ( f a b s ( z)> a c c ) ;
p r i n t f ( ” r o o t=%f , f n=%f , a c c=%f \n” , x3 , f ( x3 ) , a c c ) ;
break ;
}}}
float f ( float x)
{ f l o a t a , b , s , h , k ; i n t i ,N=800; a = 0 . 0 ; b=p i ;
h=(b−a ) /N;
s =0.0;
f o r ( i =1; i<=N−1; i +=2)
{ s+=4∗g ( x , a+i ∗h)+2∗g ( x , a+( i +1)∗h ) ; }
s=s+g ( x , a)+g ( x , b ) ; s ∗=(h / ( 3 ∗ 2 ∗ p i ) ) ;
return ( s ) ; }

The lowest root obtained is at 2.723019.

5. Use the integral representation of the Bessel function:

Page 76 of 77;
INTEGRATION COMPUTER PROGRAMMING-2015

Z2π
1
Jn (z) = cos(z cos(x))dx

0

to find its zeroes in the range 0 ≤ z ≤ 12 by secant method.

#include <s t d i o . h>


#include <math . h>
#define g ( x , z ) ( c o s ( x∗ c o s ( z ) ) )
#define p i 3 . 1 4 1 5 9
main ( )
{
f l o a t x , xmin=0 ,xmax =12.0 , x i n c =0.5 , t e s t , x l , x1 , x2 , x3 , a c c =0.00001 , z ;
i n t f a c t ( i n t n ) ; f l o a t f ( f l o a t x ) ; i n t l s =0;
f o r ( x=xmin ; x<=xmax ; x+=x i n c )
{ t e s t=f ( x ) ∗ f ( x+x i n c ) ;
i f ( t e s t <0)
{ l s =1; x1=x ; x2=x+x i n c ;
do
{ x3 = ( x1 ∗ f ( x2)−x2 ∗ f ( x1 ) ) / ( f ( x2)− f ( x1 ) ) ;
x1 = x2 ; x2 = x3 ; z = f a b s ( x1−x2 ) ; l s ++;
}
while ( f a b s ( z)> a c c ) ;
p r i n t f ( ” r o o t=%f , f n=%f , a c c=%f \n” , x3 , f ( x3 ) , a c c ) ;
}}}
float f ( float x)
{ f l o a t a , b , s , h , k ; i n t i ,N=800; a = 0 . 0 ; b=2∗ p i ;
h=(b−a ) /N;
f o r ( k = 1 . 0 ; k <=12.0; k+=0.2)
{ s =0.0;
f o r ( i =1; i<=N−1; i +=2)
{ s+=4∗g ( x , a+i ∗h)+2∗g ( x , a+( i +1)∗h ) ; }
s=s+g ( x , a)+g ( x , b ) ; s ∗=(h / ( 3 ∗ 2 ∗ p i ) ) ;
} return ( s ) ; }

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;

You might also like