MATLAB

Return

Introduction

Help

To access help files you can type help, help functionname, or use the help menu.

Saving Files

There are many commands, which you can type into the Matlab Command Window, but more complicated routines are best saved to an ASCII or text file as filename.m. Then you can run filename, or call up a user-defined function, more readily. Most editing can be done in Notepad. Be sure to save as a a file with *.m and not *.txt or *.m.txt.

Your files should be saved in a directory, which is in the path Matlab is using. You can see the paths by typing path in the Command Window. You can add a new path to the list. This you will probably need to do often when working in the students labs. For example, when working off a floppy, you should first enter

path(path,'a:\')

For a different directory, replace 'a:\' with the appropriate directory structure.

1D Maps

One can do simple 1D maps, like the logistic map.

Create a file logistic.m in the bin directory and type logistic in the Command Window.

clear
N=400;
x(1)=0.4;
r0=3.8;
for n=2:N,
    x(n)=r0*x(n-1)*(1-x(n-1));
end
plot(x,'b')

Here is a template using a function in Matlab.

Create a function file map1d.m

function y=map1d(x,r)
y=r*x*(1-x);

Create the iteration file

N=100;
x(1)=0.4;
r=2.1;
for n=2:N
x(n)=map1d(x(n-1),r);
end

Things one can do with the output:

  • List orbits: x
  • List single element: x(100)
  • change format: format long
  • Plot orbit: plot(x)
  • Format plot: Title, axis, points, colors, .... These things you can read about in the tutorials and the help pages.

 

2D Maps

One can now think about 2D maps, like the henon map. This is not the most efficient use of Matlab, but shows the basic apsects of programming.

Create a file called henon.m

N=500;
a=0.5;
b=.3;
axis([-2.5 2.5 -2.5 2.5]);
[x0 y0 bb]=ginput(1);           (Line for mouse input)
hold
while bb<2,                           (Left click to run, right click to end)
plot(x0,y0);
for i=1:N,
    x1=a-x0^2+b*y0;
    y1=x0;
    x0=x1;
    y0=y1;
    plot(x0,y0);
    [x0,y0]
end
[x0 y0,bb]=ginput(1);
end
hold;

Another version of henon can do the iteration of the map for a set of initial conditions and plot the eigenvectors on the plot as well.

Create a file called henon2.m

a=1.4;
b=.3;
r=0.1;
clf;
axis([-2.5 2.5 -2 2]);
title(['Henon Map ' num2str(N) ' Iterations']);
hold on
for i=0:200
    x0 = .883896+r*cos(i*2*pi/200);
    y0 = .883896+r*sin(i*2*pi/200);
    plot(x0,y0);
    for i=1:N,
        x1=a-x0^2+b*y0;
        y1=x0;
        x0=x1;
        y0=y1;
        if i==N
           plot(x0,y0);
        end
    end
end
line([.883896+0,.883896+.155946],[.883896+0,.883896+1]);
line([.883896+0,.883896-1.9237],[.883896+0,.883896+1]);
hold off;

Template for 2D maps:

We can do a general 2D map like the 1D map above. Here is a map called the standard map. Put it into a fucntion file called map2d.m

function [x1,y1]=map2d(x0,y0)
k=0.7;
y1=y0-k/(2*pi)*sin(2*pi*x0);
x1=y1+x0;
y1=mod(y1,1);
x1=mod(x1,1);

We then write a routine to iterate this map and plot the points as it iterates. Save this under a name like iter2d.m and type iter2d in the Command Window.

x(1)=0.0;
y(1)=0.3;
axis([0 1 0 1]);
hold on
p = plot(x(1),y(1),'.','EraseMode','none','MarkerSize',2); 
for n=2:1000
   [x(n),y(n)]=map2d(x(n-1),y(n-1));
   set(p,'XData',x(n),'YData',y(n))
   drawnow
end
hold off
A modification of this routine will allow us to plot orbits for various initial conditions. We add ginput and a while loop.
axis([0 1 0 1]);
[x(1) y(1) bb]=ginput(1);
hold on
p = plot(x(1),y(1),'.','EraseMode','none','MarkerSize',2); 
while bb<2,
   for n=2:500
      [x(n),y(n)]=map2d(x(n-1),y(n-1));
      set(p,'XData',x(n),'YData',y(n))
      drawnow
   end
   [x(1) y(1) bb]=ginput(1);
end
hold off

Iterated Function Systems

Iterated functions systems are just another example of iteration. As an example, here is the IFS code for generating a fern and watching the evolution of the plot. Again, create a file with a name like fern.m. 

y=[0,0];
p = plot(y(1),y(2),'.','EraseMode','none','MarkerSize',2); 
N=5000;
xa=-5;
xb=5;
ya=0;
yb=12;
axis([xa xb ya yb]);
hold on
a=[0,.85,0.2,-0.15];
b=[0,.04,-.26,.28];
c=[0,-.04,.23,.26];
d=[.16,.85,.22,.24];
e=[0,0,0,0];
f=[0,1.6,1.6,.44];
Pr=[.01,.85,.07,.07];
x0 = 0;
y0 = 0;
for i=1:N
    r = rand;
    if r < Pr(1) 
       x1 = a(1) * x0 + b(1) * y0 + e(1);
       y1 = c(1) * x0 + d(1) * y0 + f(1);
    elseif  r < Pr(1) + Pr(2) 
       x1 = a(2) * x0 + b(2) * y0 + e(2);
       y1 = c(2) * x0 + d(2) * y0 + f(2);
    elseif r < Pr(1) + Pr(2)+Pr(3) 
       x1 = a(3) * x0 + b(3) * y0 + e(3);
       y1 = c(3) * x0 + d(3) * y0 + f(3);
    else
       x1 = a(4) * x0 + b(4) * y0 + e(4);
       y1 = c(4) * x0 + d(4) * y0 + f(4);
    end 
    x0 = x1;
    y0 = y1;
    set(p,'XData',x0,'YData',y0)
    drawnow
 end
hold off

A slightly more compact form of this code is given below. Notice the use of matrices and commented lines.

x=[0;0];
p = plot(x(1),x(2),'.','EraseMode','none','MarkerSize',2); 
N=5000;
axis([-5 5 0 12]);
hold on
% Enter each affine transformation Ax+T
A1=[0 0; 0 .16];
T1=[0;0];
A2=[.85 .04; -.04 .85];
T2=[0;1.6];
A3=[0.2 -.26; .23 .22];
T3=[0;1.6];
A4=[-.15 .28; .26 .24];
T4=[0;.44];
% Enter the probabilities
Pr=[.01,.85,.07,.07];
% The iteration
for i=1:N
    r = rand;
    if r < Pr(1) 
       y = A1*x+T1;    
    elseif  r < Pr(1) + Pr(2) 
       y = A2*x+T2; 
    elseif r < Pr(1) + Pr(2)+Pr(3) 
       y = A3*x+T3; 
    else
       y = A4*x+T4; 
    end 
    x=y;
    set(p,'XData',y(1),'YData',y(2))
    drawnow
 end
hold off

Systems of ODEs

We will see shortly that systems of ODEs describe dynasmical systems and have all of the wonderful types of orbits we have been studying. You might want to numerically solve such systems. Here is an example of solving a system of two first order differential equations. The system is

                  x'(t) = y(t)
                  y'(t) = - x(t).

First create a function file, like odefunc.m as shown below.

function dy = odefunc(t,y)
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = -y(1);

Now you can use one of many ode solvers provided by Matlab, or write your own. We will use ode45, Type on the Command line

[T,Y] = ode45('odefunc',[0 1000],[2 0]);

To get a phase portrait, type plot(Y(:,1),Y(:,2)), or to get a solution vs time, type plot(T,Y(:,1)).

New! - ODEs

To plot a family of solutions of a first order ODE, one can use the following set of files. The first file, which can be saved as ode1d.m, calls the differential equation from the function file odefun1.m. Clicking on the graph produces a solution curve for the solution passing through the point (t0,y(t0). 

clf
axis([0 100 -1 2]);
[t0 x0 bb]=ginput(1);           
hold on
while bb<2, 
   [T,Y] = ode45('odefun1d',[t0 100],x0); 
   plot(T,Y), 
   [t0 x0,bb]=ginput(1);
end
hold off;
function dy = odefun1d(t,y)
dy = 0;    
dy = y*(1-y);

New! - ODEs

Here is a routine for generating a variety of initial condtions for a system of differential equations specified in a file odefunc.m. Both the phase portrait and the soltion are plotted in two windows. In this example a system needs to be input. An example of a linear systems follows.

% Routine to plot phase portrait and solutions to system odefunc
% Plot window h1 is the solution window
% Plot window h2 is the portrait window
% Note the new functions for plotting windows:
%      subplot, title
% 
h1=subplot('Position',[0.1 0.03 0.8 0.18]);
axis([0 100 -1 1]);          % axis for solution
title(['Solution vs Time']);
subplot(h1)
hold on
h2=subplot('Position',[0.15 0.27  0.5 0.7]);
axis([-5 5 -5 5]);           % axis for portrait
title(['y(2) vs y(1)']);
subplot(h2)
[x0 y0 bb]=ginput(1);           
hold on
while bb<2, 
   [T,Y] = ode45('odefunc',[0 100],[x0 y0]); 
   subplot(h2)
   plot(Y(:,1),Y(:,2)),       % plot portrait
   subplot(h1)
   plot(T,Y(:,1)),            % plot solution
   subplot(h2)
   [x0 y0,bb]=ginput(1);
end
hold off;
subplot(h1)
hold off
Here is the set of function for a linear system of ODEs.
function dy = odefunc(t,y)
dy = zeros(2,1);    % a column vector
a=-1.0;
b=-1.0;
c=1.0;
d=0.5;
dy(1) = a*y(1)+b*y(2);
dy(2) = c*y(1)+d*y(2);

Analysis of Times Series

In the following code the sunspot data, which comes with Matlab, is plotted using a delay plot. Saving the program as sundata.m, it is a function, which can be called where the argument m is the delay step; e.g, sundata(2) will plot the data x(n) vs x(n-m) for n>m-1.

function [x,y]=sundata(m)
load sunspot.dat
year=sunspot(:,1);
data=sunspot(:,2);
[Nrow Ncol]=size(data);
for i=1:Nrow-m
    j=i+m;
    x(i)=data(j,1);
    y(i)=data(i,1);
end    
plot(x,y)
In the next ocde, one can enter a set of data and perform a delay plot. 
function embed1(data)
[Nrow Ncol]=size(data);
for i=1:Ncol-1
    j=i+1;
    x(i)=data(1,j);
    y(i)=data(1,i);
 end    
 plot(x,y)
One can also generate their own periodic data:
function [x,y]=genperdat(N,amp,per,cf)
for i=1:N
   y=rand(1,N)*cf;
end   
x=linspace(0,20,N);
[Nrow Ncol]=size(per);
for i=1:Ncol
   y=amp(i)*cos(2*pi/per(i)*x)+y;
end

% Type of command code
% [xz,yz]=gendat(500,[1 -0.5],[0.5 0.6],0); plot(xz,yz)

Mandelbrot Set

Here is code to produce the Mandelbrot set. This is mostly for showing how to compute the set in a naive way. It is not as fast as other codes, like Fractint. 

function mand(a,b,g,N)
a0=real(a);
b0=imag(a);
a1=real(b);
b1=imag(b);
axis([a0 a1 b0 b1]);
hold on
p = plot(a0,b0,'.','EraseMode','none','MarkerSize',2); 
p = plot(a1,b1,'.','EraseMode','none','MarkerSize',2); 
dr=(a1-a0)/g(1);
di=(b1-b0)/g(2);
for i1=1:g(1)
   xr0= a0+dr*i1;
   for j=1:g(2)
      x0= (b0+di*j)*i+xr0;    
      x=0;
      fl=0;
      for i=1:N
         x=x^2+x0;
         if abs(x)>2
            fl=i;
         	break   
         end
      end
      if fl==0 
          set(p,'XData',real(x0),'YData',imag(x0))
          drawnow
      end
   end  
end  
hold off

% Sample Command to type
% mand(-2-2i,2+2i,[500 500],100)
Typing mand(-2-2i,2+2i,[500 500],100) generated this figure.