A problem using central divided difference error order

This is a problem I asked in the first examination of my Numerical Methods course in Spring 2009.   The question is that if one gives you an approximate value of the derivative of a function at a certain point using the central divided difference formula for two different step sizes, would you be able to find a better estimate of the derivative?

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://www.youtube.com/numericalmethodsguy.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Audiovisual Lectures on Numerical Methods

We have started to put audiovisual lectures on the web for an undergraduate course in numerical methods. These videos are available at the youtube site http://www.youtube.com/numericalmethodsguy.

Some of the videos are longer than 10 minutes, so they could not be loaded to the youtube site. Right now, these videos and the links to the youtube videos are available at http://nm.mathforcollege.com/videos/index.html. We are editing the longer than 10-minute videos so that they can be uploaded to youtube.

The above webpage is a better webpage to find the videos of your choice as the videos are given in a sequential fashion for a particular topic.

By end of Spring 2009, we should have uploaded all the audiovisual lectures for a typical numerical methods course.

Drop me a comment on what you think about this new avenue to learn numerical methods.

_____________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://www.youtube.com/numericalmethodsguy.  

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Proper modeling needs to precede numerical solutions

Given below is a problem where it is shown that proper modeling is key to solving engineering problems.  An improper model and a proper model may lead to opposite conclusions


_________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com and the textbook on Numerical Methods with Applications available from the lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Numerical Methods Book Printed

The blog is restarting after a month of hiatus due to winter break at my institution. As part of the dissemination plan of our NSF funded grant, we have printed a numerical methods book. This book consists of most of the numerical methods taught in a typical numerical methods course.

The book is being used in Spring 2009 in the numerical methods course at the University of South Florida and at the Arizona State University.

You can buy the textbook: Numerical Methods with Applications at the lulu storefront.

A dedicated website for the book has also been developed.

___________________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com and the textbook on Numerical Methods with Applications available from the lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Is a square matrix strictly diagonally dominant?

A square matrix A is strictly diagonally dominant if for all rows the absolute value of the diagonal element in a row is strictly greater than than the sum of absolute value of the rest of the elements in that row.

In this posting, I show a MATLAB program that finds whether a square matrix is strictly diagonally dominant by using two different methods. These are academic ways to reinforce programming skills in a student.

The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor). The html file showing the mfile and the command window output is also available.

%% IS A GIVEN SQUARE MATRIX STRICTLY DIAGONALLY DOMINANT?
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 25, 2008
% Abstract: This program shows you two ways of finding out
% if a square matrix is diagonally dominant. A square matrix is
% diagonally dominant if for all rows the absolute value of the
% diagonal element in a row is strictly greater than than the sum
% of absolute value of the rest of the elements in that row
clc
clear all
disp(‘This program shows you two ways of finding out’)
disp(‘if a square matrix is diagonally dominant. A square matrix is’)
disp(‘diagonally dominant if for all rows the absolute value of the’)
disp(‘diagonal element in a row is strictly greater than than the sum’)
disp(‘of absolute value of the rest of the elements in that row’)
disp(‘ ‘)
%% INPUTS
% The square matrix
A=[-12 1 -7 2;1 3.4 1.1 1.1; 1 0 -4.5 0;10 1 1 10];
disp (‘INPUTS’)
disp(‘Here is the square matrix’)
A
disp(‘ ‘)

%% FIRST SOLUTION
% This is based on finding for how many rows the condition
% the absolute value of the diagonal element in a row is
% strictly greater than than the sum of absolute value
% of the rest of the elements in that row.

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% count = for how many rows is the inequality met that
% the absolute value of the diagonal element in a row is
% strictly greater than than the sum of absolute value
% of the rest of the elements in that row
count=0;
for i=1:1:n
sumrow=0;
for j=1:1:n
if i~=j
sumrow=sumrow+abs(A(i,j));
end
end
if abs(A(i,i))>sumrow
count=count+1;
end
end
disp(‘FIRST WAY’)
if count==n
disp(‘Matrix is strictly diagonal dominant’)
else
disp(‘Matrix is NOT strictly diagonal dominant’)
end

%% SECOND SOLUTION
% This is based on finding for if for any row the condition
% the absolute value of the diagonal element in a row is
% strictly greater than than the sum of absolute value
% of the rest of the elements in that row is NOT met

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% flag = keeps track if the condition is not met
% flag = 1 if matrix is strictly diagonally dominant
% flag = 2 if matrix is not strictly diagonally dominant

% Assuming matrix is strictly diagonally dominant
flag=1;
for i=1:1:n
sumrow=0;
for j=1:1:n
if i~=j
sumrow=sumrow+abs(A(i,j));
end
end
% As soon as the condition is not met, it is not a strictly
% diagonally dominant matrix
if abs(A(i,i))<=sumrow
flag=2;
break;
end
end
disp(‘ ‘)
disp(‘SECOND WAY’)
if flag==1
disp(‘Matrix is strictly diagonal dominant’)
else
disp(‘Matrix is NOT strictly diagonal dominant’)
end

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com and the textbook on Numerical Methods with Applications available from the lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Skipping numbers in picking the lotto numbers

In a previous post, a reader wanted a program that would skip some numbers that he thinks would presumably not show up because they were picked in the previous week’s lotto.  Although the probability of the same numbers being picked is the same, human psychology and the lack of true randomness is a viable factor in writing such a program.  So here it is where you can input what numbers you do not want to see picked up by the program.

The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is also available.

%% PICKING THE LOTTO NUMBERS WHILE SKIPPING SOME
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 10, 2008
% Abstract: This program chooses randomly m unique numbers to play
% the lotto.  Lotto numbers allowed are positive integers from
% xlow to xhigh, while numbers asked to be skipped are not chosen
clc
clear all
disp(‘This program chooses randomly m unique numbers to play’)
disp(‘the lotto.  Lotto numbers allowed are positive integers from’)
disp(‘xlow to xhigh, while numbers asked to be skipped are not chosen’)
disp(‘  ‘)
%% INPUTS
% xlow= lowest integer allowed
xlow=1;
% xhigh = highest integer allowed
xhigh=53;
% number of integers to be picked
m=6;
% numbers to be skipped
lottoskip=[ 2  5  7 12  21];
disp (‘INPUTS’)
fprintf(‘Lowest integer allowed=%g’,xlow)
fprintf(‘\nHighest integer allowed=%g’,xhigh)
fprintf(‘\nNumbers to be picked=%g\n\n’,m)
disp(‘Lotto numbers skipped’)
lottoskip
disp(‘ ‘)
% Using the random number generator rand to get the lotto numbers.
% rand generates numbers between 0 and 1.  So we multiply that by xhigh-xlow+1
% and shift it by xlow and then floor it to get the integer part.
%% SOLUTION
% number of lotto numbers to be skipped from being chosen
mskip=length(lottoskip);
i=1;
while (i<=m)
lottonum(i)=floor(xlow+rand*(xhigh-xlow+1));
% flag= variable that keeps track of previous numbers
% being same or different from the number picked
% also it checks if the numbers asked to skipped are
% not chosen
flag=0;
% checking against previous numbers chosen
for j=1:1:i-1
if (lottonum(i)==lottonum(j))
flag=1;
end
end
% checking against numbers to be skipped
for j=1:1:mskip
if (lottonum(i)==lottoskip(j))
flag=1;
end
end
if flag==0
i=i+1;
end
end
%% OUTPUT
disp(‘ ‘)
disp(‘OUTPUT’)
disp(‘The lotto numbers picked are’)
fprintf(‘%g ‘,lottonum)
disp (‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Is a square matrix diagonal or not?

A square matrix A is diagonal if all the elements on the off-diagonal are zero. That is, A(i,j)=0 for i~=j.

In this posting, I show a MATLAB program that finds whether a square matrix is diagonal by using three different methods. These are academic ways to reinforce programming skills in a student.

The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor). The html file showing the mfile and the command window output is also available.

%% IS A GIVEN SQUARE MATRIX A DIAGONAL MATRIX?
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 15, 2008
% Abstract: This program shows you three ways of finding out
% if a square matrix is a diagonal matrix. A square matrix is
% diagonal if all the off-diagonal elements are zero, that is
% A(i,j)=0 for i~=j.
clc
clear all
disp(‘This program shows you three ways of finding out’)
disp(‘if a square matrix is a diagonal matrix.’)
disp(‘A square matrix is diagonal if all the off-diagonal’)
disp(‘elements are zero, that is A(i,j)=0 for i~=j.’)
disp(‘ ‘)
%% INPUTS
% The square matrix
A=[1 0 0 0;0 3.4 0 0; 0 0 -4.5 0;0 0 0 0];
disp (‘INPUTS’)
disp(‘Here is the square matrix’)
A
disp(‘ ‘)

%% FIRST SOLUTION
% This is based on counting the number of zeros on
% off the diagonal. If this count is n^2-n then it
% is a diagonal matrix, otherwise it is not a diagonal matrix

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% count = how many zeros not on the diagonal
count=0;
for i=1:1:n
for j=1:1:n
if A(i,j)==0 & i~=j
count=count+1;
end
end
end
disp(‘FIRST WAY’)
if count==n^2-n
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

%% SECOND SOLUTION
% This is based on finding if any of the off-diagonal elements
% are nozero. As soon as this condition is met, the matrix can be
% deemed not diagonal. If the condition is never met, the matrix is
% diagonal

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);
% flag = keeps track if it is diagonal or not
% flag = 1 if matrix is diagonal
% flag = 2 if matrix is not diagonal

% Assuming matrix is diagonal
flag=1;
for i=1:1:n
for j=1:1:n
% flag=2 if off-diagonal element is nonzero.
if A(i,j)~=0 & i~=j
flag=2;
end
end
end
disp(‘ ‘)
disp(‘SECOND WAY’)
if flag==1
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

%% THIRD SOLUTION
% This is based on finding if the sum of the absolute value of
% the off-diagonal elements is nonzero.
% If the sum is nonzero, the matrix is NOT diagonal.
% If the sum is zero, the matrix is diagonal

%size gives how many rows and columns in the A matrix
rowcol=size(A);
n=rowcol(1);

% sum_off_diagonal= sum of absolute value of off-diagonal elements
sum_off_diagonal=0;
for i=1:1:n
for j=1:1:n
if i~=j
sum_off_diagonal=sum_off_diagonal+abs(A(i,j));
end
end
end

disp(‘ ‘)
disp(‘THIRD WAY’)
if sum_off_diagonal==0
disp(‘Matrix is diagonal’)
else
disp(‘Matrix is NOT diagonal’)
end

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Picking lotto numbers

In these tough economic times, more people are playing lotto to hit it big.  We have written a MATLAB program that chooses unique and random numbers for a lotto.  The inputs are the number of numbers to pick, and the range in which the numbers can be picked.

The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is also available.

Can you modify this program to put the numbers in an ascending order?

How does the random number generator work?

How can I get random numbers that are different every time I run the program?

%% PICKING THE LOTTO NUMBERS
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 10, 2008
% Abstract: This program chooses randomly m unique numbers to play
% the lotto.  Lotto numbers allowed are positive integers from
% xlow to xhigh
clc
clear all
disp(‘This program chooses randomly m unique numbers to play’)
disp(‘the lotto.  Lotto numbers allowed are positive integers from’)
disp(‘xlow to xhigh’)
disp(‘  ‘)
%% INPUTS
% xlow= lowest integer allowed
xlow=1;
% xhigh = highest integer allowed
xhigh=53;
% number of integers to be picked
m=6;

disp (‘INPUTS’)
fprintf(‘Lowest integer allowed=%g’,xlow)
fprintf(‘\nHighest integer allowed=%g’,xhigh)
fprintf(‘\nNumbers to be picked=%g’,m)
disp(‘ ‘)
% Using the random number generator rand to get the lotto numbers.
% rand generates numbers between 0 and 1.  So we multiply that by xhigh-xlow+1
% and shift it by xlow and then floor it to get the integer part.
%% SOLUTION
i=1;
while (i<=m)
lottonum(i)=floor(xlow+rand*(xhigh-xlow+1));
% flag= variable that keeps track of previous numbers
% being same or different from the number picked
flag=0;
for j=1:1:i-1
if (lottonum(i)==lottonum(j))
flag=1;
end
end
if flag==0
i=i+1;
end
end
%% OUTPUT
disp(‘ ‘)
disp(‘OUTPUT’)
disp(‘The lotto numbers picked are’)
fprintf(‘%g ‘,lottonum)
disp (‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Comparing two series to calculate pi

Many series are used to calculate the value of pi.  In this blog, we compare two series, one by Gregory and another by Ramanujan.

Here is a MATLAB program that does the comparison for you.  The MATLAB program can be downloaded as a Mfile (better to download it, as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is also available.

%% COMPARING TWO SERIES FOR VALUE OF PI
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : October 30, 2008
% Abstract: This program compares results for the value of
% pi using a) Gregory series and b) Ramanajun series
clc
clear all
clf
format long
disp(‘This program compares results for the value of’)
disp(‘pi using a) Gregory series and b) Ramanajun series’)
disp(‘  ‘)
disp(‘Gregory series’)
disp(‘pi=sum over k from 0 to inf of (4*((-1)^k/(2*k+1))’)
disp(‘  ‘)
disp(‘Ramanajun Series’)
disp(‘1/pi=sum over k from 0 to infinity of 2*sqrt(2)/9801*((4k)!*(1103+26390k)/(k!)^4*396^(4*k))’)

%% INPUTS.
%If you want to experiment this the only parameter
% you should and can change.
% Maximum number of terms
n=30;

%% PROGRAM

%% GREGORY SERIES
pi_gregory=0;
for i=1:1:n
pi_gregory=pi_gregory+(-1)^(i+1)*4*(1/(2*i-1));
pi_gregory_array(i)=pi_gregory;
end

%% RAMANUJAN SERIES
pi_ram=0;
for i=0:1:n-1
pi_ram=pi_ram+2*sqrt(2)/9801.0*(factorial(4*i))*(1103.0+26390.0*i)/((factorial(i)^4)*(396)^(4*i));
pi_ram_array(i+1)=1/pi_ram;
end

%% THE OUTPUT
disp(‘ ‘)
fprintf(‘\nNumber of Terms = %g’,n)
fprintf(‘\nGregory Series Value = %g’,pi_gregory)
fprintf(‘\nRamanujan Series Value = %g’,1/pi_ram)
disp( ‘   ‘)

%% PLOTTING THE TWO SERIES AS A FUNCTION OF TERMS
x=1:1:n;
hold on
xlabel(‘Number of terms’)
ylabel(‘Value of pi’)
title(‘Comparing Gregory and Ramanujan series’)
plot(x,pi_gregory_array,’color’,’blue’,’LineWidth’,2)
hold on
plot(x,pi_ram_array,’color’,’black’,’LineWidth’,2)
legend(‘Gregory Series’,’Ramanajun Series’,1)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

An automatic integrator using Trapezoidal rule

How would you know how many segments to use in a Trapezoidal rule of integration to get an accurate value of the integral?  This can be done by applying the Trapezoidal rule for 1 segment rule, then 2 segment rule, followed by 4 segment rule and so on.  As soon as the absolute relative approximate error (page 5-6) between the consecutive answers becomes less than the pre-specified tolerance chosen by the user, you would have your integral within the accuracy you desired.

Here is a MATLAB program that does that for you.  The MATLAB program that can be downloaded at http://nm.mathforcollege.com/blog/trapezoidal_rule_automatic.m (better to download it as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is here: http://nm.mathforcollege.com/blog/html/trapezoidal_rule_automatic.html

% Simulation : Using Trapezoidal rule as an automatic integrator

% Language : Matlab 2007a

% Authors : Autar Kaw, http://nm.mathforcollege.com

% Mfile available at
% http://nm.mathforcollege.com/blog/trapezoidal_rule_automatic.m

% Last Revised : October 12, 2008

% Abstract: This program uses multiple-segment Trapezoidal
% rule to integrate f(x) from x=a to x=b within a pre-specified tolerance

clc
clear all

disp(‘This program uses multiple-segment Trapezoidal rule as an automatic integrator’)
disp(‘to integrate f(x) from x=a to x=b within a pre-specified tolerance’)
disp(‘ ‘)
disp(‘Author: Autar K Kaw.’)
disp(‘http://autarkaw.wordpress.com’)
disp(‘http://nm.mathforcollege.com’)
disp(‘ ‘)

%INPUTS.  If you want to experiment, these are the only variables
% you should and can change.
% a = Lower limit of integration
% b = Upper limit of integration
% nmax = Maximum number of segments
% tolerance = pre-specified tolerance in percentage
% f = inline function as integrand
a=5.3;
b=10.7;
nmax=20000;
tolerance=0.005;
f=inline(‘exp(x)*sin(2*x)’);

% SIMULATION
disp(‘INPUTS’)
func=[‘     The integrand is =’ char(f)];
disp(func)
fprintf(‘     Lower limit of integration, a= %g’,a)
fprintf(‘\n     Upper limit of integration, b= %g’,b)
fprintf(‘\n     Maximum number of segments, nmax = %g’,nmax)
fprintf(‘\n     Pre-specified percentage tolerance, eps = %g’,tolerance)
disp(‘  ‘)
disp(‘  ‘)

% Doing the automatic integration
% Calculating the integral using 1-segment rule
previous_integral=(b-a)/2*(f(a)+f(b));
% Initializing ea as greater than pre-specified tolerance for loop to work
ea=2*tolerance;
% Starting with 2-segments inside the while loop
n=2;
while (ea>tolerance) & (n<=nmax)
h=(b-a)/n;
% Keeping track of used number of segments
nused=n;
current_integral=0;
for i=1:1:n-1
current_integral=current_integral+f(a+i*h);
end
current_integral=2*current_integral+f(a)+f(b);
current_integral=(b-a)/(2*n)*current_integral;
% Calculating the absolute relative approximate error
ea = abs((current_integral-previous_integral)/current_integral)*100;
previous_integral=current_integral;
% Doubling the number of segments for next estimate of the integral
n=n*2;
end

disp(‘OUTPUTS’)
fprintf(‘      Number of segments used  =%g’, nused)
fprintf(‘\n      Approximate value of integral is =%g’,current_integral)
fprintf(‘\n      Absolute percentage relative approximate error =%g’, ea)
if (ea>tolerance)
disp(‘  ‘)
disp(‘  ‘)
disp(‘     NOTE: The value of integral is not within the pre-specified tolerance’)
end
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.