Matlab Tutorial

MATLAB tutorial 1: using the command line

Start matlab: open a terminal and type “matlab &”

At the matlab command prompt type a = 1

and press enter. Next b = 1.0

and press enter. To see what variable are currently defined type whos

You should see something like:

>> whos

 

Name

Size

Bytes Class Attributes

a

1x1

8 double

b

1x1

8 double

An array is a list of numbers assigned to one variable. Type x = [ 0 1 2 3 4 5 ]

we may also use y = 0:0.5:5

Try

x(3) = 22

To prevent matlab from printing out the entire array every time we change something, we add a semicolon at the end of the command:

x(4) = -10;

Arrays can be manipulated using matrix algebra: z = 2*y

or

w = 2*y+z Now try

x*x

x.*x

x*x'

x'*x

Note that the single dash after an array means “transpose” – i.e. a swap of the rows and columns. Note also that the dot in front of the multiplication sign of the second command means we multiply element-by-element.

Let's try to plot something: plot(x,3*x)

A new plot window should have appeared.

Now try

           u = x.*x

           plot(x,u)

We can save this plot as a png file by print  -dpng

or as an encapsulated postscript file by print  -depsc2

These commands created the files “figure1.png” and “figure1.eps”, which we could rename (e.g. “mv figure1.png plot_x-squared.png” on the linux command line) if we wanted to keep them.

We will now read in some data from the file “partdiv.txt”. Looks what's inside this file by typing “less partudiv.txt” on the LINUX command line. Use the up and down arrows and page-up & page-down to navigate the file. When you're done press “q”.

Load the data into matlab with:

            [congress, year, total, dems, reps] = textread('partydiv.txt', '%s %f %f %f %f',78,'headerlines',3)

Be careful here, there is a space before each % sign.

Let's make some plots:

           plot(year,seats)

           xlabel('Year')

           ylabel('Total House Seats')

Now try

          plot(year,dems./total,'-b') 

          xlabel('Year')

          ylabel('Percentage')

We can overplot by turning the "hold on"

          hold on

          plot(year,reps./total,'-r')

          hold off

We can add a legend with

          legend('democrats', 'republicans')

You can move the legend around using the GUI by first clicking the arrow, then selecting the legend box and dragging it.

You can get a sense for how correlated (or anti-correlated) a dataset is by looking at the “correlation coefficient”:

corrcoef(dems./total,reps./total)

gives you a 2x2 matrix of coefficients relating each column to itself (along the diagonal) and to the other column(s) (off-diagonals).

TASK

Now that you have completed your tutorial, let’s do some science!

Inside your home directory, you should find a file called “RecentIndices.txt”. Use Linux to see what is inside the file and what every column of the data is. The first column is the year and the second is the month of the observation. The third column is the observed sunspot number, the 8th is the observed radio flux, and the 10th is the observed AP index (a geomagnetic index).

To complete this task, you must

1)generate an image (png) files that contain the sunspot number as a function of time, the radio flux as a function of time, and the observed AP index as a function of time. Plot the lines in different colors and make a legend

2)generate 2 postscript files, one that contains the scatter plot of the sunspot number on the x-axis and radio flux on the y axis and one that contains the sunspot number vs AP index. Calculate the correlation coefficient between the two values and include that information in the title of the plot.

Be sure to label all axes of your plots.

When you are done with this task, email the results (3 files) to jacob.heerikhuisen@uah.edu

MATLAB tutorial 2: writing a script

Start matlab: open a terminal and type “matlab &”

In terminal type: “emacs first.m &” to open a text editor

In the text editor type:

clear; % clears all variables so each run starts fresh x = 0:0.1:10;

y = x.^2; plot(x,y)

save this file and run it by typing “first” on the matlab command line. This should open a plot in a new window.

Some points to remember from the first tutorial:

x is an “array” that takes on the values {0, 0.1, 0.2, 0.3, ….., 9.8, 9.9, 10.0}

the semicolon at the end of each line prevents matlab from printing out the result (try to remove one of the semicolons).

We had to put a dot in front of the power symbol to tell matlab to multiply x by itself “element-by-element\;, rather than by using matrix algebra.

Now change the second line of your code to y = x^n;

save and run it.

This should have given you an error that “n” in undefined. Fix this by adding a line n = 3

above the place where n is first used. Save and run it again.

You can also turn the program into a function, with “n” as an input variable. Try adding function [] = first(n)

to the top of the program. Save and run “first” from the command line. Now try “first(3)” from the command line.

Let's write a function to plot some 2D data. First we will create a separate function that takes in x and y values and generates a 2D grid of z values:

function z = func_2D(x,y)

z = sin(2*pi*x) .* cos(3*pi*y);

save this file as “func_2D.m”. Then create another file “second.m” function [] = second(n,m)

dx = 2/n; dy = 2/m; x = -1:dx:1; y = -1:dy:1;

grid_x = meshgrid(x,y); grid_y = -meshgrid(y,x)';

z = func_2D(grid_x,grid_y); pcolor(x,y,z)

pause

surf(x,y,z)

shading interp

Loops are useful for manipulating a list (array) of data one element at a time. Some of the above tricks of using matrix multiplication are more efficient, but loops are useful for complicated data. Consider the program:

n = 10;

data = rand(n,1); heads = 0;

tails = 0; for k = 1:n

dat = data(k); if (dat > 0.5)

heads = heads + 1;

else

tails = tails + 1;

end

end

disp(['fraction of heads = ', num2str(heads/n)]) disp(['fraction of tails = ', num2str(tails/n)])

TASK: Write a program to generate a histogram

A histogram is a way of sorting and describing date. For instance, if there are 10 people in the room and you sort them into ages (3 20 year olds, 4 21 year olds, and 3 22 year olds) you have made a histogram. The input to the program would be the ages of people, i.e. ages= [20, 21, 22, 21, 20, 20, 21, 22, 22, 21], the output would be the array that describes the ages (bin= [20,21,22]) and the numbers of people in each age bin (hist =[3,4,3]). In your program, you should have at least 1 input (like ages) and at least 2 outputs (like bin and hist). The other thing you have to consider is the size of the bins. For instance, what if your bin was in 10 year increments, bin = [10, 20, 30]? All the ages would be in the 20 bin (h = [0,10,0]) and your histogram would be pretty boring. Alternatively, the bin size could be too small, you could have increments of days or even minutes. In that case, the bin array would be huge and there would only be one person in each one bin so h would be mostly 0 with 10 1s. So binsize could also be an input into your program or you should figure out a way to make the program to calculate the correct binsize.

Once you have written you histogram program, try it out on some real data. Make a histogram of the sunspot number data we looked at in the first tutorial. Make a histograms of the number of months as a function of sunspot number.

Write the plot into a eps file and send it to jacob.heerikhuisen@uah.edu along with your histogram program.

MATLAB tutorial 3: writing codes to solve a problem

An Orbit Around the Sun

We will try to compute the orbit of an object (assumed massless) around the Sun. Open a new text file and start by clearing all variables:

clear;

Next let's define the length of an astronomical unit (AU) and the number of seconds in a day:

au = 1.49598e11

% m

day = 60*60*24

% s

We are going to integrate the equations of motion in thexy-plane: dxdt =vdydt =dvdt=advdt=ay

The acceleration can be computed from the gravitational force:

GM r

r̂

where the vector r is a unit vector, and the mass of the object cancels out. A unit vector is a vector of length 1, obtained by taking the position vector and dividing it by its length |r|. Using this information, and the fact that =( x , ygives:

a=(a

, a

)=

GM (x , y )

x

y

r2

r

or

ax=GM x ;a y=

GM y

r3

r3

The four differential equations require four initial conditions in order to have a unique solution. Let's choose the following:

xx = 1 * au; yy = 0; vx = 0; vy = 29800;

As we integrate over the orbit, we need to keep track of the position of the object, so that we can plot out the orbit once it's been calculated. It's most efficient to “declare” the arrays where we will store the x & y positions up-front:

npoints = 1000;

x_vect = zeros(1,npoints); y_vect = zeros(1,npoints);

Finally, we need to choose a time-step for the integration procedure. dt = 1*day

We are going to use the simplest integration method, known as “Euler's method”. The basic idea is that at each point in time you take a small step in the direction of the local velocity, and then update the velocity according to the local force. The rest of the program looks something like:

x_vect(1) = xx/au; y_vect(1) = yy/au; for kk = 2:npoints

[ax,ay]=orbit_func(xx,yy,vx,vy);

xx = xx + vx*dt;

yy = yy + vy*dt;

vx = vx + ax*dt;

vy = vy + ay*dt;

x_vect(kk) = xx/au; y_vect(kk) = yy/au;

endplot(x_vect,y_vect,'-k')axis([-1.5 1.5 -1.51.5]) axis square

Here we make use of an external function that computes the accelerations at a given point: function [ax,ay]=orbit_func(xx,yy)

G = 6.6742e-11;

% m^3 kg^{-1} s^{-2}

M_s = 1.9889e30;

% kg

r = sqrt(xx^2+yy^2);

F_factor = -G*M_s / r^3;

ax = F_factor*xx;

ay = F_factor*yy;

Run the program to see the orbit. Try changing the step size. Is a smaller step size more accurate? How can you tell?

The above method can be greatly improved by using a more accurate method. The next simplest is known as “second order Runge-Kutta”. The method involves taking the average of the current velocity with an estimate of the “future” velocity:

xx1 = xx + vx*dt; vx1 = vx + ax*dt;

xx = xx + 0.5*(vx + vx1)*dt;

Similar terms are required for the y position and the x & y velocities – the latter requiring additional accelerations ax1 & ay1 at the positions xx1 & yy1.

Compare the accuracy of the Euler & RK2 methods for a number of orbits and a number oftime-step sizes.

TASKS

Use your code to answer the following:

1.Suppose a comet has a perihelion of 0.2 AU, and an aphelion of 100AU. What speed will the comet have at perihelion?

2.Suppose we want to send a spacecraft to Mars (assume a distance of 1.5 AU from the Sun). What is the minimum launch speed relative to the earth that could get the spacecraft to Mars.

3.What is the escape speed for the solar system?

Try the next two tasks if you have time.

4.Extend your code to two finite masses orbiting around each other.

5.Make a movie of the two orbiting masses. The following code will be useful to create frames for your movie:

outname=['frame'num2str(time)'.png'] 

print(gcf,'-dpng',outname)