MATLAB programming – Page 2 – Numerical Methods Guy

Differentiating a Discrete Function with Equidistant Points

Many students ask me how do I do this or that in MATLAB. This is a new addition to the “How do I do that in MATLAB” series.

In this blog, I show you how to find the first derivative of a discrete function y(x). We are assuming that the x values are eqidistant and the data is sorted in ascending or descending order by the x values. The latter requirement can be relaxed easily in programs such as MATLAB where one can use the sortrows command to put the data in the required order.

To keep the accuracy of all calculated first derivatives to be the same, we use the following formulas:

For the first data point, we use the forward divided difference formula

f ‘(x) =(-f(x+2h) + 4 f(x+h) -3 f(x))/(2h)+order(h^2)

For the interior points, we use the central divided difference formula

f ‘(x) =(f(x+h) -f(x-h))/(2h)+order(h^2)

For the last data point, we use the backward divided difference formula

f ‘(x) =(f(x-2h) – 4 f(x-h) +3 f(x))/(2h)+order(h^2)

Here are the links for the program.

The mfile is here

The published version of the mfile 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 find the first derivative of a discrete function y(x) if the
% x values are equidistant.
%% SUMMARY
% Language : Matlab 2010a;
% Authors : Autar Kaw and Sri Garapati;
% Mfile available at
% http://nm.mathforcollege.com/blog/discrete_diff_equidistant_blog.m
% Last Revised : January 17, 2012;
% Abstract: This program shows you how to differentiate discrete data if
% the x values are equally spaced
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘ This program shows you how to differentiate discrete data if’)
disp(‘ the x values are equally spaced ‘)

disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar Kaw and Sri Garapati of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/discrete_diff_equidistant_blog.m’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   January 17, 2012’)
disp(‘ ‘)

%% INPUTS

% Inputs assuming
%    that three or more points are given
%    all x values are equidistant
%    x values are in ascending or descending order.
x=[2.3   3.4   4.5    5.6   6.7   7.8];
y=[4.6   7.9   13.0   12.3 3.2   1.9];
%% DISPLAYING INPUTS
disp(‘ ‘)
disp(‘INPUTS’)
% Creating a matrix to print input data as a table
data=[x;y]’;
disp(‘   X Data     Y Data’)
% Printing the input data as a table
disp(data)

%% THE CODE

% n returns the number of data points
n=length(x);
% delta is the distance between consecutive x values
delta=x(2)-x(1);

% “dy” is an array which stores the value of the derivative at each x-value

% finding the derivative from the 2nd order accurate forward divided
% difference formula at the first data point
dy(1)=(-y(3)+4*y(2)-3*y(1))/(2*delta);

% finding the derivative from the 2nd order accurate central divided
% difference formula at the second to (n-1)th data point
for i=2:1:n-1
dy(i)=(y(i+1)-y(i-1))/(2*delta);
end

% finding the derivative from the 2nd order accurate backward divided
% difference formula at the first data point
dy(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*delta);

% creating a matrix with input data and calculated derivative value at
% each data point for printing as a table
xdy=[x’ y’ dy’];

%% DISPLAYING OUTPUTS
disp(‘ ‘)
disp(‘OUTPUTS’)
disp(‘ XData YData Derivative’)
% printing the input data points and calculated derivative values(Outputs)
disp(xdy)
_________________________________________________

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, the textbook on Introduction to Programming Concepts Using MATLAB, and the YouTube video lectures available at http://nm.mathforcollege.com/videos. Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Reading an excel file in MATLAB

Recently I taught a volunteer class to professional engineers on MATLAB. Two of the most requested items of interest were
1. How do I read an excel file?
2. How do I do curve fitting?

We address the first question here. It is easy to read an excel file with the xlsread command but what do you with it once the file has been assigned. So we took a simple example of an excel spreadsheet where the first column consists of a student number and the second column has the examination scores of the students. You are asked to find the highest score.

The MATLAB program link is here.
The HTML version of the MATLAB program is here.
The Excel file used in the MATLAB program is here

It is better to download (right click and save target) the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or you can read the html version for clarity and sample output.

%% READING AN EXCEL SPREADSHEET IN MATLAB
% Language : Matlab 2008a
% Authors : Autar Kaw
% Last Revised : December 12, 2010
% Abstract: This program shows you how to read an excel file in MATLAB
% The example has student numbers in first column and their score in the
% second column
clc
clear all
disp(‘This program shows how to read an excel file in MATLAB’)
disp(‘Matlab 2008a’)
disp(‘Authors : Autar Kaw’)
disp(‘Last Revised : December 12, 2010’)
disp(‘http://nm.mathforcollege.com’)
disp(‘ ‘)

%% INPUTS
% We have two column data and it has headers in the first row.
% That is why we read the data from A2 to B32.
A=xlsread(‘c:\users\grades.xls’,’A2:B32′);
disp (‘The data read from the excel spreadsheet is’)
disp(A)
disp(‘ ‘)
%% SOLUTION
% Finding the number of rows and columns
sizem=size(A);
rows_A=sizem(1);
cols_A=sizem(2);
% Assigning the scores to a vector called score
for i=1:1:rows_A
    score(i)=A(i,2);
end
% Using the max command to find the maximum score
% HW: Write your own function “max”
maxscore=max(score);
% Finding which student got the highest score
for i=1:1:rows_A
   if score(i)==maxscore
       student_no=i;
       break;
   end
% HW: What if more than one student scored the highest grade??
end
%% OUTPUT
disp(‘ ‘)
disp (‘OUTPUT’)
fprintf(‘Student Number# %g scored the maximum score of %g’,…
    student_no,maxscore)
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,
the textbook on Introduction to Programming Concepts Using MATLAB, and
the YouTube video lectures available at http://nm.mathforcollege.com/videos

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

Using int and solve to find inverse error function in MATLAB

In the previous post, http://autarkaw.wordpress.com/2010/08/24/finding-the-inverse-error-function/, we set up the nonlinear equation to find the inverse of error function. Using the int and solve MATLAB commands, we write our own program to find the inverse error function.

It is better to download (right click and save target) the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or you can read the html version for clarity and sample output.

%% FINDING INVERSE ERROR FUNCTION
% In a previous blog at autarkaw.wordpress.com (August 24, 2010),
% we set up a nonlinear equation to find the inverse error function.
% In this blog, we will solve this equation.
% The problem is given at
% http://nm.mathforcollege.com/blog/inverseerror.pdf
% and we are solving Exercise 1 of the pdf file.

%% TOPIC
% Finding inverse error function

%% SUMMARY

% Language : Matlab 2010a;
% Authors : Autar Kaw;
% Mfile available at
% http://nm.mathforcollege.com/blog/inverse_erf_matlab.m;
% Last Revised : August 27, 2010
% Abstract: This program shows you how to find the inverse error function
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to’)
disp(‘   find the inverse error function’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/inverse_erf_matlab.m’)
disp(‘ ‘)
disp(‘PROBLEM STATEMENT’)
disp(‘   http://nm.mathforcollege.com/blog/inverseerror.pdf’)
disp(‘        Exercise 1’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   August 27, 2010’)
disp(‘ ‘)

%% INPUTS
% Value of error function
erfx=0.5;

%% DISPLAYING INPUTS

disp(‘INPUTS’)
fprintf(‘ The value of error function= %g’,erfx)
disp(‘ ‘)
disp(‘ ‘)

%% CODE
syms t x
inverse_erf=solve(int(2/sqrt(pi)*exp(-t^2),t,0,x)-erfx);
inverse_erf=double(inverse_erf);
%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘ Value of inverse error function from mfile is= %g’,inverse_erf)
fprintf(‘\n Value of inverse error function using erfinv is= %g’,erfinv(erfx))
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,
the textbook on Introduction to Programming Concepts Using MATLAB, and
the YouTube video lectures available at http://nm.mathforcollege.com/videos

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

Solving a polynomial equation for the longest mast problem?

In the previous post, http://autarkaw.wordpress.com/2010/06/10/a-real-life-example-of-having-to-solve-a-nonlinear-equation-numerically/, we set up the polynomial equation for the problem of finding the length of the mast before which it buckles under its own weight. In this blog, we show you how the polynomial equation is solved. Since the polynomial equation has infinite terms, we also show you how we choose how many terms of the polynomial need to be used.

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

%% FINDING THE SMALLEST POSITIVE ROOT OF A POLYNOMIAL EQUATION FOR A
% In a previous blog at autarkaw.wordpress.com (June 10), we set up a
% polynomial equation that would allow us to find the longest mast
% that can be setup before it buckles under its own weight.
% In this blog, we will find the root of this equation.
% The problem is given at
% http://nm.mathforcollege.com/blog/length_of_mast.pdf
% and we are solving Exercise 1 of the pdf file.

%% TOPIC
% Smallest positive root of a polynomial equations with infinite terms

%% SUMMARY

% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://nm.mathforcollege.com/blog/rootinfinite.m;
% Last Revised : July 4, 2010
% Abstract: This program shows you how to find the smallest positive
% real root of a polynomial equations with infinite terms
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to’)
disp(‘   find the smallest positive real root’)
disp(‘   of a polynomial equations with infinite terms’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://nm.mathforcollege.com/blog/rootsinfinite.m’)
disp(‘ ‘)
disp(‘PROBLEM STATEMENT’)
disp(‘   http://nm.mathforcollege.com/blog/length_of_mast.pdf)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   July 4, 2010’)
disp(‘ ‘)

%% INPUTS
% prespecified tolerance, eps
eps=0.000001;
% maximum number of terms of polynomial
nmax=100;

%% DISPLAYING INPUTS

disp(‘INPUTS’)
fprintf(‘ The max number of terms of the polynomial chosen, nmax= %g’,nmax)
fprintf(‘\n The prespecified tolerance, eps= %g’,eps)
disp(‘ ‘)
disp(‘ ‘)

%% CODE
for N=2:1:nmax
    % The above looop is to see how many terms we should take of the
    % infinite polynomial
aa(1)=-3.0/8.0;
    % Setting up the polynomial via recursive relations
for i=2:1:N
aa(i)=-3*aa(i-1)/(4*i*(3*i-1));
end
% Since it is a polynomial of order N,
% there are N+1 coefficients
% To set up the polynomial for MATLAB
% N+1 th coefficient is the constant term
% N th coefficient is the term of order 1
% and so on till 1st coefficient is of order N
bb(N+1)=1;
for i=1:1:N
    bb(N-i+1)=aa(i);
end

% Finding all the roots of the Nth order polynomial
abc=roots(bb);

% Finding the first real positive root so that it
% would be used as the starting minimum value available
for i=1:1:N
    if isreal(abc(i))==true & abc(i)>0
        minval=abc(i);
        break;
    end
end

% Finding the smallest positive real root

for i=1:1:N
    if isreal(abc(i))==true & abc(i)>0
        if (abc(i) < minval)
            minval=abc(i);
        end
    end
end
% Checking if prespecified tolerance is met
if N>2
    absea=abs((minval-previous)/minval)*100;
    if absea<=eps
        terms_needed=N;
        break;
    end
end
previous=minval;
end

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘ The number of terms used in the polynomial is %g’,terms_needed)
fprintf(‘\n The smallest positive real root is %g’,minval)
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://www.youtube.com/numericalmethodsguy.

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

A short online quiz for the MATLAB conditional statements

Frequent testing has been proven to be effective in learning. Here we have an online quiz on MATLAB conditional statements (only if-then-else), where some of the questions are calculated (meaning that the numbers in these questions change when you retake the quiz). So give it a try.
http://nm.mathforcollege.com/EML3041/studymate/MATLABifelseend.htm

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Transforming 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.

A short online quiz on the for-end loops in MATLAB

Frequent testing has been proven to be effective in learning. Here we have an online quiz on MATLAB loops (just the for-end loops in this quiz), where some of the questions are calculated (meaning that the numbers in these questions change when you retake the quiz). So give it a try and soon we will be adding questions on other programming basics.
http://nm.mathforcollege.com/EML3041/studymate/MATLABforendloops.htm

____________________________________________________

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.

A short online quiz on MATLAB basics

Frequent testing has been proven to be effective in learning. Here we have an online quiz on MATLAB basics, where some of the questions are calculated (meaning that the numbers in these questions change when you retake the quiz). So give it a try and soon we will be adding questions on programming basics.
http://nm.mathforcollege.com/EML3041/studymate/matlabbasics.htm

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.

MATLAB code for bubble sort

In the previous blog, we spelled out the bubble sort algorithm for putting an array of numbers in an ascending order. In this post, I am posting the matlab program. It is better to download the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or see the html version for clarity and sample output.

%% PUTTING AN VECTOR OF NUMBERS IN AN ASCENDING ORDER?
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : November 8, 2009
% Abstract: This program shows you how to put a vector
% of numbers in an ascending order using the bubble sort method
clc
clear all
disp(‘This program shows the bubble sort method’)
disp(‘to put a vector of numbers in an ‘)
disp(‘ascending order’)
disp(‘Matlab 2007a’)
disp(‘Authors : Autar Kaw’)
disp(‘Last Revised : November 8, 2009’)
disp(‘http://nm.mathforcollege.com’)
disp(‘ ‘)
%% INPUTS
% The vector of numbers
disp (‘INPUTS’)
disp(‘Input the vector of numbers’)
A=[18 7 6 15 4 13];
disp(A)
%% SOLUTION
% Number of entries, n
n=length(A);
% making (n-1) passes
for j=1:1:n-1
    % comparing each number with the next and swapping
    for i=1:1:n-1
    if A(i)>A(i+1);
        % temp is a variable where the numbers are kept
        % temporarily for the switch
        temp=A(i);
        A(i)=A(i+1);
        A(i+1)=temp;
    end
    end
end

%% OUTPUT
disp(‘ ‘)
disp (‘OUTPUT’)
disp (‘The ascending matrix is’)
disp(A)

_______________________________________________________

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.

Bubble sorting

Last week, I was teaching how to randomly pick lotto numbers using MATLAB. The problem was that some of the numbers that were getting picked were identical. We solved this by using comparisons until the current number picked is different from the previously selected numbers. We blogged on this a few months ago. But there is still an aesthetic problem of how the numbers are presented. The numbers are not in an ascending or descending order. This is a good time to show how to do this using the simplest (if not the most efficient) procedure called the bubble sort.

Let’s suppose someone asks you to put [8 7 9 5 4] in an ascending order.

Starting from the first number, you compare the number with the next number, and see if it is greater. If it is, you swap the numbers. You continue to do this with the second number, third number and so on until the second last number. What this does is bubble the largest number to the end.

[8 7 9 5 4 ] -> [7 8 9 5 4] (as 8>7) -> [7 8 9 5 4] (as 8 is not > 9)-> [7 8 5 9 4] (as 9 is >5) -> [7 8 5 4 9] (as 9>4). See how the largest number is at the end.

Now repeat this. [7 8 5 4 9] -> [7 8 5 4 9] (as 7 is not >8] -> [7 5 8 4 9] (as 8>5) -> [7 5 4 8 9] (as 8>4) -> [7 5 4 8 9 ] (as 8 is not >9)

Now repeat this. [ 7 5 4 8 9] -> [5 7 4 8 9] -> [ 5 4 7 8 9] -> [ 5 4 7 8 9] -> [ 5 4 7 8 9]

Now repeat this. [5 4 7 8 9] -> [ 4 5 7 8 9] -> [ 4 5 7 8 9] -> [ 4 5 7 8 9] -> [ 4 5 7 8 9]

It looks like we are done. If n is the number of the numbers in the array, it takes (n-1) swaps within each of the (n-1) repetitions. So we do not have to guess how many swaps it takes or how many repetitions it takes.

To make the bubble sort efficient, we can do the following: 1) Since with each repetition the largest number bubbles up, we may need to do less swaps. For the first repetition, we will do (n-1) swaps, for the next repetition, we do the first (n-2) swaps, and so on. 2) We can also keep track of number of swaps taking place in a repetition. If no swaps take place in a repetition, no more repetitions are needed.

_______________________________________________________

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.

How do I numerically solve an ODE in MATLAB?

The other day a student came to ask me for help in solving a second order ordinary differential equation using the ode45 routine of MATLAB. To use ode45, one needs to be familiar with how the inputs are required by MATLAB. The understanding of these inputs is important to use ode45 successfully in problems that are more complex than solving a second order ODE.

The ordinary differential equation was
2y”+3y’+5y=7 exp(-x), y(0)=11, dy/dx(0)=13
This has to put in the state variable form by reducing it by using
y’=z
That gives
y’=z with the corresponding initial conditions as y(0)=11
Then
2y”+3y’+5y=7 exp(-x)
reduces to
2z’ + 3z+5y=7exp(-x)
z’ =(7exp(-x)-3z-5y)/2 with the corresponding initial conditions as z(0)=13

So as needed by MATLAB, call y as y(1) and z as y(2)
dy(1)=y(2), y(1) at x=0 is 11
dy(2)=(7exp(-x)-3y(2)-5y(1))/2, y(2) at x=0 is 13

These equations are now put in a MATLAB function we call odestate.m
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=(7*exp(-x)-3*y(2)-5*y(1))/2;

To solve the ODE, the
The inputs are
1) the function odestate
2) The outputs are required between x=0 and x=17,
hence entered as [0 17]
3) The initial conditions are y(0)=11 and dy/dx(0)=13,
hence entered as [11 13]

The outputs are
1) X= array of x values between 0 and 17
2) Y= matrix of 2 columns;
first column is the y(x)
second column is dy/dx(x)
The MATLAB code then is
[X,Y]=ode45(@odestate,[0 17],[11 13]);

Click the links for the MATLAB mfiles for the function odestate.m and the ODE solver odetest.m

_________________________________________________________________________________

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