matlab – Page 2 – Numerical Methods Guy

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.

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

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

How do I read data from a textfile 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 read data from an external data file. We will use the simplest command to read from a text file, that is, the fgets command.

The MATLAB program link is here and the text file is here.

So let us look at an example. Let’s suppose someone asks you to find the dot product of two vectors: A=[7 11 13 23] and B=[3 5 17 29]. The program for that simply is

A=[7 11 13 23]
B=[3 5 17 29]
n=length(A)
dot_product=0
for i=1:1:n
dot_product=dot_product+A(i)*B(i)
end

What if now the data for the vectors are given in an external text file? How does MATLAB read the data and assign it properly to A and B.

The first thing is to make the text file. Let’s save this data in the file called /blog/vectordata.txt. The text file is here.
7 11 13 23
3 5 17 29

Now how is this data read from the data file.
First you tell MATLAB what the name of the file is
filen=’W:/blog/vectordata.txt’

Now the file needs to be opened by using the fopen command. This also assigns an integer to the file called the filehandler. This integer is automatically assigned by MATLAB and is a unique integer for any open file. The ‘r’ stands for the file being opened for reading. The other choices are ‘w’ for writing and ‘a’ for appending.
fh=fopen(filen,’r’)

Now we need to get the data from the file which is identified through the file identifier variable. The fgets command does this and every time fgets is used, it reads the next line.

A_vector=fgets(fh)

However, this reads the whole line of data as a single string. How do we now get it to be a vector A with the 4 entries. We do this by using the sscanf command that parses the string variable into the format used by sscanf. The format used here is %f which is a floating point format. The other formats for reading numbers include %e for scientific format and %d for integers.
A=sscanf(A_vector,’%f’)

Similarly, one can now read the B vector via fgets command and parse it using the sscanf command.

The MATLAB program link is here and the text file is here.

__________________________________________

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)

_______________________________________________________

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.

_______________________________________________________

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

_________________________________________________________________________________

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

The continue statement in MATLAB

The continue statement in MATLAB is used to pass control to the next iteration in for and while statements. Let’s suppose someone wants to find and print the value of k^2-50 for all integers in [-10,10] domain. The mfile for that is given below.

% For integers k=-10,-9,….,9,10,
% the function k^2-50 will take positive as
% well as negative values.
%For example, for k=-9, k^2-50=31; for k=1,
% k^2-50=-49; for k=8, k^2-50=14.
% The loop below will calculate values of k^2-50 for
% all values of requested k.
for k=-10:1:10
val=k^2-50;
fprintf(‘\n k=%g val=%g’,k,val)
end

___________________________________________________________

Let’s suppose now you are asked to calculate and print value of k^2-50 for all integers in [-10,10] domain but only if (k^2-50) is positive.

% The loop below will calculate and print values of k^2-50
% for all values of the requested k
% for which k^2-50 is positive.
for k=-10:1:10
    if (k^2-50<0)
        continue;
    end
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
end

____________________________________________________________

Can you do what you did above using the while statement. Yes, the MATLAB code is given below.

% Equivalent in while
% The loop below will calculate values of k^2-50
% for all values of the requested k for which k^2-50 is positive.
k=-10;
while (k<=10)
    if (k^2-50>0)
        val=k^2-50;
        fprintf(‘\n k=%g val=%g’,k,val)
    end
k=k+1;
end

_______________________________________________________

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

The break statement in MATLAB

The break statement in MATLAB is used to break out of a loop – a for or while statement, that is, it terminates the execution of the loop. Let’s suppose someone wants to find the value of k^2-50 for all integers in [-10,10] domain. The mfile for that is given below.

% For integers k=-10,-9,….,9,10,
% the function k^2-50 will take positive as
% well as negative values.
%For example, for k=-9, k^2-50=31; for k=1,
% k^2-50=-49; for k=8, k^2-50=14.
% The loop below will calculate values of k^2-50 for all values of requested k.
for k=-10:1:10
val=k^2-50;
end

___________________________________________________________

Let’s suppose now you are asked to calculate value of k^2-50 for all integers in [-10,10] domain but only until k^2-50 becomes negative.

% The loop below will calculate values of k^2-50
% for all values of the requested k
% until it turns negative
for k=-10:1:10
    if (k^2-50<0)
        break;
    end
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
end

____________________________________________________________

Can you do what you did above using the while statement. Yes, the MATLAB code is given below.

% Equivalent in while
% The loop below will calculate values of k^2-50
% for all values of the requested k until it turns negative
k=-10;
while (k<=10) & (k^2-50>0)
    val=k^2-50;
    fprintf(‘\n k=%g val=%g’,k,val)
    k=k+1;
end

_______________________________________________________

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