fishy population

% fishy.m
%
% usage:
% fishy
%
% purpose:
% Load some data on fish populations
%
% data provided:
% fish: fish(:,1) age in years
% fish(:,2) mean weight (mass) in gms
% fish(:,3) current no./sq. km of a given age
% fish(:,4) no. of eggs (/1000) /per female of given age
%

% Some data

% ages
fish(:,1)=[1 2 3 4 5 6 7 8 9 10]';
% weights
fish(:,2)=[120 250 350 500 680 820 980 1110 1220 1300]';
% frequencies/sq km
fish(:,3)=0.5*[10000, 8197, 2725, 907, 302, 101, 33, 11, 4, 1]';
% egg production
fish(:,4)=20*[0, 0.7, 1.4, 2.2, 3.0, 3.8, 4.4, 5.0, 5.0, 5.0]';

disp(' fish data loaded in matrix fish(:,:)')
disp(' ')

M-file for producing a random permutation of 1,...,n and then breaking it into disjoint cycles

% M-file for producing a random permutation of 1,...,n and then breaking it into disjoint cycles

n=input('Type n ');

p=randperm(n)

pstart=p; % Just in case you want to be reminded what p was at the end

disp('Now follow the cycles of this permutation')

i=0;

while i i=i+1;
if p(i)>0
v=[];
flag=0; % flag is just to trick
% the computer into entering
% the while loop
first=i;
while (first ~=i) | (flag==0)

% ~= means not equal
% and | means or
flag=1;
r=first;
v=[v r];
first=p(first);
p(r)=0; % The entries in the cycle are changed
% to zero to indicate that they
% have been `used'
end;
disp(v)
end;
end;

"cobweb" diagram for Mobius sequences

% Plots a "cobweb" diagram for Mobius sequences, starting value x0. You have to input the constants a,b,c,d and the lower and upper values of x and y which will appear on the graphics screen.
%
% Iteration is done 30 times, or until it seems that convergence is definitely taking place.
%
% If you break into the program with , then remember to type "hold off" to release the graph. You can still see the graphics window.


a=input('Type the value of a ');
b=input('Type the value of b ');
c=input('Type the value of c ');
d=input('Type the value of d ');
x0=input('Type the value of x0 ');

xl=input('Type the lower limit of x ');
xu=input('Type the upper limit of x ');
yl=input('Type the lower limit of y ');
yu=input('Type the upper limit of y ');


dx=(xu-xl)/100;
x=xl:dx:xu;
y=(a*x+b)./(c*x+d);
plot(x,y)
axis([xl xu yl yu]);

hold on

plot([xl xu],[xl xu],'b') % Draws the diagonal in blue

clear x y iterate
x(1)=x0;
y(1)=0;
x(2)=x0;
y(2)=0;
x(3)=x0;
y(3)=(a*x0+b)/(c*x0+d);

iterate(1)=x0;
iterate(2)=y(3);

plot(x,y,'g')
i=1;
x(1)=x0+.1; % To trick it into entering this while loop
while abs(x(3)-x(1)) > .001 & i < 30
x(1)=x(3);
y(1)=y(3);
x(2)=y(1);
y(2)=x(2);
x(3)=x(2);
y(3)=(a*x(2)+b)/(c*x(2)+d);
iterate(i)=y(3);
plot(x,y,'g')
pause
i=i+1;
end

xlabel('No more iterations will be done. Press to return to MatLab')
pause

iterate'

hold off

Cubic graph

% M-file for drawing some cubic graphs y=x^3-ax+1 You need to input the lower and upper bounds of x and y on the screen, and also the value of a.

clg

xl=input('Type the lower bound of x ');
xu=input('Type the upper bound of x ');
yl=input('Type the lower bound of y ');
yu=input('Type the upper bound of y ');



a=input('Type the value of a ');

xstep=(xu-xl)/200;


x=xl:xstep:xu;
y=x.*x.*x - a*x + ones(size(x));

plot(x,y,'r')
axis([xl xu yl yu]);
hold on

u=[xl xu];
v=[0 0];
plot(u,v,'w')

hold off

Newton-Raphson iteration for a quadratic equation with complex roots p, q

% M-file to do complex Newton-Raphson iteration for a quadratic equation with complex roots p, q

i=sqrt(-1); % This is just in case i has been used as an integer lately

p=input('Type the complex number p, using the form p1 + i*p2 ');
q=input('Type the complex number q, using the form q1 + i*q2 ');

a=-p-q;
b=p*q;

% So p and q are the roots of the equation x^2+ax+b=0

hold off

%**********************************************************

xl=input('Type the lower bound of x ');
xu=input('Type the upper bound of x ');
yl=input('Type the lower bound of y ');
yu=yl+xu-xl;

%**********************************************************

% The next section draws crosses at the two roots.
% The cross size can be adjusted via the next line.

cross_size=(xu-xl)/30;

p1=real(p); p2=imag(p);
q1=real(q); q2=imag(q);


x=[p1-cross_size p1+cross_size];
y=[p2 p2];
plot(x,y,'r')

hold on

x=[p1 p1];
y=[p2-cross_size p2+cross_size];
plot(x,y,'r')

x=[q1-cross_size q1+cross_size];
y=[q2 q2];
plot(x,y,'g')

y=[q2-cross_size q2+cross_size];
x=[q1 q1];
plot(x,y,'g')

%**** Axis commands here

axis([xl,xu,yl,yu]);
axis('square');

% The screen is now a square area in the Argand diagram,
% with xl + i yl at bottom left and xu + i yu at top right.

%******************************************************

% Now we iterate using the Newton-Raphson
% method, looking for the roots p, q of z^2 + az + b = 0,
% starting at 500 random points z = x + iy.
% If the iterates approach p we mark them in red and
% if they approach q we mark them in green. If after
% 20 iterates they are still far from both p and q
% then we abandon that starting point.

epsilon=0.00001;

for j=1:500
x=xl+rand*(xu-xl);
y=yl+rand*(yu-yl);

count=0;
z=x+i*y;
while (abs(z-p)>epsilon)&(abs(z-q)>epsilon)&(count<20)
z=(z^2-b)/(2*z+a);
count=count+1;
end;
if abs(z-q)<=epsilon
plot(x,y,'go');
elseif abs(z-p)<=epsilon
plot(x,y,'ro');
end;
end;

hold off


i=sqrt(-1); % This is just in case i has been used as an integer lately

p=input('Type the complex number p, using the form p1 + i*p2 ');
q=input('Type the complex number q, using the form q1 + i*q2 ');

a=-p-q;
b=p*q;

% So p and q are the roots of the equation x^2+ax+b=0

hold off

%**********************************************************

xl=input('Type the lower bound of x ');
xu=input('Type the upper bound of x ');
yl=input('Type the lower bound of y ');
yu=yl+xu-xl;

%**********************************************************

% The next section draws crosses at the two roots.
% The cross size can be adjusted via the next line.

cross_size=(xu-xl)/30;

p1=real(p); p2=imag(p);
q1=real(q); q2=imag(q);


x=[p1-cross_size p1+cross_size];
y=[p2 p2];
plot(x,y,'r')

hold on

x=[p1 p1];
y=[p2-cross_size p2+cross_size];
plot(x,y,'r')

x=[q1-cross_size q1+cross_size];
y=[q2 q2];
plot(x,y,'g')

y=[q2-cross_size q2+cross_size];
x=[q1 q1];
plot(x,y,'g')

%**** Axis commands here

axis([xl,xu,yl,yu]);
axis('square');

% The screen is now a square area in the Argand diagram,
% with xl + i yl at bottom left and xu + i yu at top right.

%******************************************************

% Now we iterate using the Newton-Raphson
% method, looking for the roots p, q of z^2 + az + b = 0,
% starting at 500 random points z = x + iy.
% If the iterates approach p we mark them in red and
% if they approach q we mark them in green. If after
% 20 iterates they are still far from both p and q
% then we abandon that starting point.

epsilon=0.00001;

for j=1:500
x=xl+rand*(xu-xl);
y=yl+rand*(yu-yl);

count=0;
z=x+i*y;
while (abs(z-p)>epsilon)&(abs(z-q)>epsilon)&(count<20)
z=(z^2-b)/(2*z+a);
count=count+1;
end;
if abs(z-q)<=epsilon
plot(x,y,'go');
elseif abs(z-p)<=epsilon
plot(x,y,'ro');
end;
end;

hold off

Newton-Raphson iteration for a quadratic equation with complex roots p, q

% M-file to do complex Newton-Raphson iteration for a quadratic equation with complex roots p, q
% This M-file takes 10 random starting points and traces the paths
% obtained by joining up successive approximations until
% these get very close to one of the roots.

i=sqrt(-1); % This is just in case i has been used as an integer lately

p=input('Type the complex number p, using the form p1 + i*p2 ');
q=input('Type the complex number q, using the form q1 + i*q2 ');

a=-p-q;
b=p*q;

% So p and q are the roots of the equation x^2+ax+b=0


%**********************************************************

xl=input('Type the lower bound of x ');
xu=input('Type the upper bound of x ');
yl=input('Type the lower bound of y ');
yu=yl+xu-xl;

%**********************************************************

% The next section draws crosses at the two roots.
% The cross size can be adjusted via the next line.

cross_size=(xu-xl)/30;

p1=real(p); p2=imag(p);
q1=real(q); q2=imag(q);


x=[p1-cross_size p1+cross_size];
y=[p2 p2];
plot(x,y,'r')

hold on

x=[p1 p1];
y=[p2-cross_size p2+cross_size];
plot(x,y,'r')

x=[q1-cross_size q1+cross_size];
y=[q2 q2];
plot(x,y,'g')

y=[q2-cross_size q2+cross_size];
x=[q1 q1];
plot(x,y,'g')

% Axis commands come next

axis([xl,xu,yl,yu]);
axis('square');
% The screen is now a square area in the Argand diagram,
% with xl + i yl at bottom left and xu + i yu at top right.


%******************************************************

% Now we iterate using the Newton-Raphson
% method, looking for the roots p, q of z^2 + az + b = 0,
% starting at 10 random points z = x + iy.
% The computer draws a broken line joining the starting
% point to the successive approximations to a root.
% You hit Enter to choose the next random starting point.

epsilon=0.00001;

for j=1:10
x=xl+rand*(xu-xl);
y=yl+rand*(yu-yl);

count=0;
z=x+i*y;
while (abs(z-p)>epsilon)&(abs(z-q)>epsilon)&(count<20)
znew=(z^2-b)/(2*z+a);
plot([real(z) real(znew)], [imag(z) imag(znew)])
z=znew;
count=count+1;
end;
pause
end;

hold off

Chop

function Xout=chop(Xin,n)

% Usage: y = chop(x,n)

% rounds elements of x to n decimal digits.

% For example, chop(24.567,4) returns 24.57

% chop(24.564,4) returns 24.56

% chop(24.564,2) returns 25.00

% chop(24.564,1) returns 20.00

%

% Users of more recent versions of Matlab 4.1c+ should

% use the system provided command "chop" instead.



% Abort if only one input argumentis is given.

if nargin<2, help chop, return, end



% abs(Xin)=A*10^E, A=0.ddddd... & E=floor(log10(|Xin|))+1

% round the decimals in abs(Xin)/10^K & K=E-n.

zero = (Xin==0);

Xout=abs(Xin)+zero; size_x = size(Xout);

E_X= (10*ones(size_x)).^(floor(log10(Xout))-n+1);

Xout=round(Xout./E_X).*E_X; Xout=sign(Xin).*Xout.*(1-zero);

Exponential distribution illustration

% c6exp.m - Chapter 6: exponential distribution illustration

help c6exp

disp('This M-file illustrates 3 probability distributions...')


rand('seed',0); n=20000; % Reset the 'seed' value
m=1।2; s=1।2; % mean and std
x = unirand(0,2।4,n,1); % Uniform in (0,2।4)
randn('seed',0); % Reset the 'seed' value
y = normrand(m,s,n,1); % Normal N(1।2,1।2)


rand('seed',0); % Reset the 'seed' value
z = exprand(m,n,1); % Exponential



subplot(311); hist(x,90) % Plot x
axis([0 2।4 0 225]) % Cut the right half

hold on; plot([m m],[0 225],':r'); hold off

subplot(312); hist(y,20) % Plot y
hold on; plot([m m],[0 4000],':r'); hold off



subplot(313); hist(z,40)% Plot z

axis([0 10 0 6852]) % Cut the right half

hold on; plot([m m],[0 6852],':r'); hold off



%%% The following lines are for "cumsum" %%%

t=cumsum(z);

fprintf('For n = %d, t first and last = %10.4f, %10.4f\n',n,t(1),t(n));

format compact; tm=hist(z,40); sum(tm);tm=max(tm)