matlab – Page 4 – Numerical Methods Guy

How do I differentiate in MATLAB?

Many students ask me how do I do this or that in MATLAB. So I thought why not have a small series of my next few blogs do that. In this blog, I show you how to differentiate a function.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

___________________________________________

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB. Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I differentiate a function?

%% SUMMARY
% Language : Matlab 2008a
% Authors : Autar Kaw
% Mfile available at
% http://nm.mathforcollege.com/blog/differentiation.m
% Last Revised : March 21, 2009
% Abstract: This program shows you how to differentiate a given function

%% INTRODUCTION
clc
clear all
disp(‘ABSTRACT’)
disp(‘   This program shows you how to differentiate’)
disp(‘   a given function and then find its value’)
disp(‘   at a given point’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/differentiation.m’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   March 21, 2009’)
disp(‘ ‘)

%% INPUTS

% Differentiate 7 exp(3*x) once and find the value of the
% first derivative at x=0.5
% Define x as a symbol
syms x
% Defining the function to be differentiated
y=7*exp(3*x);
% Defining the point where you want to find the derivative
xx=0.5;

%% DISPLAYING INPUTS
disp(‘INPUTS’)
func=[‘ The function is to be differentiated is ‘ char(y)];
disp(func)
fprintf(‘ Value of x where you want to find the derivative, x= %g’,xx)
disp(‘ ‘)
disp(‘ ‘)

%% THE CODE
% Finding the derivative using the diff command
% Argument 1 is the function to be differentiated
% Argument 2 is the variable with respect to which the
% function is to be differentiated – the independent variable
% Argument 3 is the order of derivative
dydx=diff(y,x,1);
% subs command substitues the value of x
dydx_val=subs(dydx,x,xx);
%% DISPLAYING OUTPUTS
disp(‘OUTPUTS’)
derivative_func=[‘ The derivative of function ‘ char(y) ‘ is ‘ char(dydx)];
disp(derivative_func)
fprintf(‘ Value of dydx at x=%g is =%g’,xx,dydx_val)
disp(‘ ‘)

______________________________________________________________________

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://nm.mathforcollege.com/videos and 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.

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.

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

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

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.

Another improper integral solved using trapezoidal rule

In a previous post, I showed how Trapezoidal rule can be used to solve improper integrals. The example used in the post was an improper integral with an infinite interval of integration.

In an example in this post, we use Trapezoidal rule to solve an improper integral where the integrand becomes infinite. The integral is $latex \int_{0}^{b} 1/sqrt{x} dx $. The integrand becomes infinite at x=0. Since x=0 would be one of the points where the integrand will be sought by the multiple-segment Trapezoidal rule, we choose the value of the integrand at x=0 to be zero (any other value would do too – a better assumption would be f(h), where h is the segment width in the multiple-segment Trapezoidal rule).

Here is a MATLAB program that shows you the exact value of the integral and then compares it with the multiple-segment Trapezoidal rule. The convergence is slow but you can integrate improper integrals using Trapezoidal rule.

The MATLAB program that can be downloaded at http://nm.mathforcollege.com/blog/trapezoidal_improper_sqrtx.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_improper_sqrtx.html

% Simulation : Can I use Trapezoidal rule for an improper integral?

% Language : Matlab 2007a

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

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

% Last Revised : October 8, 2008

% Abstract: This program shows use of multiple segment Trapezoidal
% rule to integrate 1/sqrt(x) from x=0 to b, b>0.

clc
clear all

disp(‘This program shows the convergence of getting the value of ‘)
disp(‘an improper integral using multiple segment Trapezoidal rule’)
disp(‘Author: Autar K Kaw. autarkaw.wordpress.com’)

%INPUTS. If you want to experiment, these are the only variables
% you should and can change.
% b = Upper limit of integration
% m = Maximum number of segments is 2^m
b=9;
m=14;

% SIMULATION
fprintf(‘\nFinding the integral of 1/sqrt(x) with limits of integration as x=0 to x=%g’,b)

% EXACT VALUE OF INTEGRAL
% integrand 1/sqrt(x)
syms x
f=1/sqrt(x);
a=0;
valexact=double(int(f,x,a,b));
fprintf(‘\n\nExact value of integral = %f’,valexact)
disp( ‘ ‘)

f=inline(‘1/sqrt(x)’);
%finding value of the integral using 16,…2^m segments
for k=4:1:m
n=2^k;
h=(b-a)/n;
sum=0;
for i=1:1:n-1
sum=sum+f(a+i*h);
end
% See below how f(a) is not added as f(a)=infinity. Instead we
% use a value of f(a)=0. How can we do that? Because as per integral calculus,
% using a different value of the function at one point or
% at finite number of points does not change the value of the
% integral.
sum=2*sum+0+f(b);
sum=(b-a)/(2*n)*sum;
fprintf(‘\nApproximate value of integral =%f with %g segments’,sum,n)
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.

The BMI (Body Mass Index) Program

In 1998, the federal government developed the body mass index (BMI) to determine ideal weights. Body mass index is calculated as 703 times the weight in pounds divided by the square of the height in inches, the obtained number is then rounded to the nearest whole number (Hint: 23.5 will be rounded to 24; 23.1 will be rounded to 23; 23.52 will be rounded to 24)

BMI Categories:

BMI <19 (underweight)
19≤BMI≤25 (healthy)
25<BMI≤30 (overweight)
BMI is >30 (obese)

Here is a MATLAB program that calculates the BMI of a person, classifies him/her in the BMI category, and then suggests a healthy weight.

It is better to download the program as single quotes in the pasted version do not translate correctly when pasted into a mfile editor of MATLAB. See the html version for clarity and sample output.

% So you want my phone number and my BMI? How shallow can you get?

% Worksheet : Finding the BMI of a person, classifying the BMI and
% suggesting a healthy weight
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : September 27, 2008

% Abstract: Finding the body mass index of a person, classifying their health
% and recommending a target weight.

%
% In 1998, the federal government developed the body mass index (BMI)
% to determine ideal weights. Body mass index is calculated as
% 703 times the weight in pounds divided by the square of the height in inches,
% the obtained number is then rounded to the nearest whole number
% (23.5 will be rounded to 24; 23.1 will be rounded to 23; 23.52 will be rounded to 24).
%
% Write a MATLAB program to do the following:
%
% Assign a value to weight in lbs, and height in inches and
% then calculate BMI as a rounded integer.
% Output a variable called health_id as
%     0 if the person’s BMI <19 (underweight)
%     1 if the person’s BMI is 19≤BMI≤25 (healthy)
%     2 if the person’s BMI is 25<BMI≤30 (overweight)
%     3 if the person’s BMI is >30 (obese)
% Output also a variable hw for healthy weight in rounded integer lbs for all the conditions.
% Use fprintf command with explanation for the inputs and outputs
clc
clear all
% INPUTS
% Weight in lbs
weight=180;
% Height in inches
height=69;

%REST OF THE PROGRAM
bmi=weight/height^2*703.0;
bmi=round(bmi);

% Assigning the proper health_id and finding the suggested healthy weight
% health_id= BMI category
% hw = healthy weight
% You can use 4 separate if-end statements or switch case statement also.
if (bmi<19)
health_id=0;
hw=19*height^2/703;
elseif (bmi>=19 & bmi<=25);
health_id=1;
elseif (bmi>25 & bmi<=30);
health_id=2;
hw=25*height^2/703;
else
health_id=3;
hw=25*height^2/703;
end
hw=round(hw);

% Printing the outputs
disp(‘The BMI MATLAB program’)
disp(‘Author: Autar K Kaw’)
disp(‘http://nm.mathforcollege.com’)
disp(‘September 28, 2008’)
disp(‘ ‘)
fprintf(‘\nWeight of person =%6.0f lbs’,weight)
fprintf(‘\nHeight of person =%6.0f inches’,height)
fprintf(‘\nYour BMI is =%g ‘,bmi)
if health_id==0
fprintf(‘\nYou are underweight. \nYour target weight is =%6.0f lbs’,hw)
elseif health_id==1
fprintf(‘\nYou are a healthy weight. \nYour target weight is =%6.0f lbs’,hw)
elseif health_id==2
fprintf(‘\nYou are overweight. \nYour target weight is =%6.0f lbs’,hw)
elseif health_id==3
fprintf(‘\nYou are obese. \nYour target weight is =%6.0f lbs’,hw)
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.

A legend used in the movie “The Happening”

Well M. Night Shyamalan may have made another disappointing movie – The Happening, but I somewhat liked it. I would give it a grade of B.

In the movie, John Leguzomo’s character, a math teacher, is distracting his fellow panicking passenger in the Jeep with a mathematical question. The question he asks her is if he gave her a penny on Day 1 of the month, two pennies on Day 2 of the month, four pennies on Day 3 of the month, and so on, how much would money would she have after a month. She shouts $300 or some odd number like that. But, do you know that the amount is actually more than a 10 million dollars (Thanks to a student who mentioned that it was a penny that John offered on the first day, not a dollar – sometimes I do feel generous).

This question is based on a story from India and it goes as follows.

King Shriham of India wanted to reward his grand minister Ben for inventing the game of chess. When asked what reward he wanted, Ben asked for 1 grain of rice on the first square of the board, 2 on the second square of the board, 4 on the third square of the board, 8 on the fourth square of the board, and so on till all the 64 squares were covered. That is, he was doubling the number of grains on each successive square of the board. Although Ben’s request looked less than modest, King Shriham quickly found that the amount of rice that Ben was asking for was humongous.

QUESTIONS:

Write a MATLAB (you can use any other programming language) program for the following using the for or while loop.

Find out how many grains of rice Ben was asking for.
If the mass of a grain of rice is 2 mg, and the world production of rice in recent years has been approximately 600,000,000 tons (1 ton=1000 kg), how many times the modern world production was Ben’s request?
Do the inverse problem – find out how many squares are covered if the the number of grains on the chess board are given to you. For example, how many squares will be covered if the number of grains on the chess board are 16?

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

Subscribe to the feed to stay updated and let the information follow you.