Finds all odd pseudoprimes to base 2 which are <= m for a value of m input by the user

% Finds all odd pseudoprimes to base 2 which are <= m for a value of m input by the user.

m=input('Enter the largest number to be tested: ');
disp('Pseudoprimes to base 2 are: ')
for n=3:2:m % Want to test odd numbers n from 3 to m.
if pow(2,n-1,n)==1 % This tests whether 2^(n-1) leaves remainder
% 1 on division by n.
% If it doesn't then we go on to the next n.
% If it does then we test whether n is in fact
% prime by trial dividing n by odd numbers
% less than the square root of n. If none of
% these divides exactly into n then n is
% certainly prime, so should not be included
% in the list of pseudoprimes.

k=3; % The first odd number to trial divide into n
while k<=sqrt(n)
if rem(n,k)==0 % Tests whether k divides exactly into n
k=n+1; % This just gets out of the while loop,
% because n has been shown to be composite
% (multiple of k) so can be displayed
% as a pseudoprime.
else k=k+2; % otherwise go on to the next odd k to trial
% divide into n.
end;
end;
if k==n+1 disp(n); end;
end; % of if pow==1 loop
end; % of for n=3:2:m loop

Power algorithm for finding rem(a^n,m)

% This implements the power algorithm for finding rem(a^n,m) (the remainder on dividing a^n by m) reliably for large numbers (up to about 8 digits).
%
% Usage: Type pow(a,n,m) where a,n,m are already defined, or are explicit numbers as in pow(2,1000,25).


function x=pow(a,n,m)
b = a;
x = 1;
while n > 0
d=rem(n,2);
if d==1
x=rem(x*b,m);
end
b=rem(b * b,m);
n=(n-d)/2;
end

M-file to find primes < 5000

% M-file to find primes < 5000

r=11;
p=[2,3,5,7];
i=4; %i counts the number of primes so far
while r<5000
s=3;
while s*s<=r & rem(r,s)>0
s=s+2;
end;
if rem(r,s)>0
i=i+1;
p(i)=r;
end;
r=r+2;
end;

Cycloid

% Approximate times of descent for a cycloid from (0,c) to (d,0), corresponding to theta=0, theta=theta1 respectively and the value r for the radius of the rolling circle.
% Uses quad8 for theta=.01 to theta1 and estimates the integral from 0 to .01 separately.

g=0.1; % g is "gravity". 0.1 seems to produce reasonable values for the time

global r_gl c_gl

c_gl=input('Type c ');
d=input('Type d ');

% The values of r and theta you need here are calculated as
% in the text, using the M-file cycfun.m and the command fsolve.

r_gl=input('Type r ');
theta1=input('Type theta1 ');

time=quad8('slide4fn',.01,theta1);
% This uses Matlab's numerical integrator to find the
% time from theta=.01 to theta=theta1.

% Now we estimate the time from theta=0 to theta=.01

theta=.01

%******************************************************************
% These are the formulae for x, y and their derivatives
% x1, y1 with respect to theta.

x=r_gl*(theta-sin(theta));
y=c_gl-r_gl+r_gl*cos(theta);
x1=r_gl*(1-cos(theta));
y1=-r_gl*sin(theta);

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

initial=sqrt(1+(x1/y1)^2)*2*sqrt(c_gl-y)/sqrt(2*g)

time=time+initial;

time

Function snowboat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%snowboat.m%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function snowboat
%
% Solver for 1st order equations associated with snowboat problem
%
% Solution is for Vector [u v x y]' where u=x' , v=y'
% It prompts for:
% mean coefft. of kinetic friction mu
% coefft of air resistance kr
% starting coords [x0 y0]
% starting veloc. [u0 v0]
% time interval dt
% and requests permission to complete another time step.
% The derivatives specifying the motion are given in snbtfn.m
%
% Uses: snbtfn.m, snowsl.m, fsnow1.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function snowboat
% Initialization
hold off; clg;format bank;
disp('Equations of motion for snow boat')
disp(' ')
disp('When you press RETURN you will see a contour plot and then')
disp(' an empty plan view of the area')
disp('You will then be asked for the initial coordinates/velocities')
disp('Press RETURN to end any pauses.')
pause;
next = 'y';
dt = 0.5; trace = 0; tol = 0.05;
%
a=2000; nx=20; dx=a/nx; % square area 2km x 2km
[X,Y]=meshdom(0:dx:a,0:dx:a); % set up contour plot of the hill
Z=snowsl(X,Y);
contour(Z,20)
title('Countour plot of snow slopes')
pause;
clg;
%
xp = [0,2000]; yp = xp;

plot(xp,yp,'.r')
title('Plan view of snow slopes')
xlabel('x')
ylabel('y')
hold on
% pause

% Input parameters
mu = input(' Input coefft of friction : ');
kr = input(' Input coefft of air res. : ');
m = 150.0; % mass of snowboat + one passenger (kgs)
k= [mu,kr,m];
while next ~= 'n'
initpos = input('Initial x,y position, enter as [x0 y0]: ');
initvel = input('Initial velocity, enter as [u0 v0]: ');
dt=input('Suggest a time step (try 10 secs if you have no better suggestion): ');
disp('Initial z and vertical velocity component are');
z0=snowsl(initpos(1),initpos(2)); [fx fy]=fsnow1(initpos(1),initpos(2));
vz0=fx*initvel(1)+fy*initvel(2);
disp([z0 vz0])

disp('please wait...')
y0 = [initpos,initvel]'; t0 = 0; tf = dt; step = 'y'; mti=1;
while step == 'y'
[tout,yout] = ode23k('snbtfn', k, t0, tf, y0, tol, trace);
%
xp = yout(:,1); yp=yout(:,2); % these are output x,y, coords
% xp = yout(:,3); yp=yout(:,4); % these are output veloc. components
%
np = length(xp); mtf=mti+np-1;
plot(xp,yp,'r');
% pause

% plot this segment of run
disp(' t x(t) y(t) u(t) v(t)')
disp([tf xp(np) yp(np) yout(np,3) yout(np,4)])
disp(' z(t) vz(t)');
z=snowsl(xp(np),yp(np)); [fx fy]=fsnow1(xp(np),yp(np));
vz=fx*yout(np,3)+fy*yout(np,4);
disp([z vz])

% store the height and speed as a function of distance from starting point
yp1(mti:mtf)=snowsl(xp,yp); %height
yp2(mti:mtf)=10.0*sqrt(yout(:,3).^2+yout(:,4).^2); %speed
xp=xp-initpos(1); yp=yp-initpos(2);
xpp(mti:mtf)=sqrt(xp.*xp+yp.*yp); %horiz. dist.
mti=mtf+1;

step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
next = input('Do you want another simulation? (y/n): ','s');
end

% now plot the height and speed overall as promised
hold off
clg
disp('Press ENTER for a plot of the height Z and the speed V ...');
disp(' as a function of the horiz. dist. travelled ');
pause

xp = [0,1500]; yp = [400,800];
plot(xp,yp,'.r')
title('Height Z is White, Speed V (x10+400) is Green')
xlabel('d')
ylabel('Z,V')
hold on
plot(xpp,yp1,'w')
yp2=yp2+400; % add 400 to make visible on same plot
plot(xpp,yp2,'g')

Species simulates the behaviour of an ecological system over a time interval

function species
%
% species simulates the behaviour of an ecological system over a time interval 'dt'। It is based on the classic Lotka-Volterra model.

% It prompts for starting populations 'X0' and 'Y0,' and permission to complete another time step.
%
% species.m uses: ode23k.m, specfn.m
%

% Initialization
disp('Ecological system with two competing species.')
disp(' ')
disp('When you press RETURN you will see an empty plot of #Y vs #X.')
disp('You will be asked for the model parameters and an initial population.')
disp('Press RETURN to end any pauses.')
pause
more = 'y';
dt = 0.1; trace = 1; tol = 0.1;
xp = [0 10]; yp = [0 10];
k = [0 0 0 0];
clg; hold off
plot(xp,yp,'.')
title('Autonomous system of competing species X and Y')
xlabel('#X')
ylabel('#Y')
hold on
pause

% Input parameters
kx = input(' Model parameters for dx/dt eqns, enter as [a b c] : ');
ky = input(' Model parameters for dy/dt eqns, enter as [d e f] : ');
k(1:3)=kx; k(4:6)=ky;
while more ~= 'n'
pop0 = input('Initial populations (from 0 to 24), enter as [X0 Y0]: ');
dt=input('Suggest a time step (try 0.1 if you have no better suggestion): ');
disp('please wait...')
y0 = pop0'; t0 = 0; tf = dt; step = 'y';
while step == 'y'
[tout,yout] = ode23k('specfn', k, t0, tf, y0, tol, trace);
xp = yout(:,1); yp=yout(:,2);
plot(xp,yp,'r'); pause
step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
more = input('Do you want another simulation? (y/n): ','s');
end
hold off

Load in tide time and height data

%%%%%%%%%%%%%%%%%%%%%%tides.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% load in tide time and height data
%
%
% Usage:-
% tides
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
echo off
load times.mat
load heights.mat
%
disp('Tide times and heights read in OK')
disp('Times in tm(1:nt), heights in hm(1:nt)')
disp(' nt dt tmin tmax ')
format bank;
nt=length(tm); dt=tm(2)-tm(1);
disp([nt,dt,tm(1),tm(nt)])
format short;

To show, for the `tomato' data, the effect of finding a line fit

% M-file = toms.m

%

% To show, for the `tomato' data, the effect of finding a line fit (not sorting data), then aquadratic fit (not sorting data) and finally a quadratic fit (sorting data). Clearly for a quadratic fit, the sorting is essential! For the linear fit it is not essential because the line, though drawn as a zigzag, still looks like a line.



tomato

C=polyfit(x,y,1);

Y=polyval(C,x);

plot(x,y,'*w',x,Y,'w')

xlabel('This is the linear fit, without sorting')



pause



C=polyfit(x,y,2);

Y=polyval(C,x);

plot(x,y,'*w',x,Y,'w')

xlabel('This is the quadratic fit, without sorting!')



pause



[sx i]=sort(x);

sy=y(i);

C=polyfit(sx,sy,2);

Y=polyval(C,sx);

plot(sx,sy,'*w',sx,Y,'w')

xlabel('This is the quadratic fit, WITH sorting')

Produce and plot a nice topographical surface

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%topog.m
%
% Produce and plot a nice topographical surface
%
% Then allows one find a local minimum of the surface
%
% Usage:-
% topog
%
% then follow the prompts.
%
% uses: snowls.m, snowmn.m, snowmx.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
echo off
clg;
a=2000; nx=20; dx=a/nx;
[X,Y]=meshdom(0:dx:a,0:dx:a);

Z=snowsl(X,Y);

mesh(X,Y,Z)

hold on
title('Snow slopes')
next=input(' More? (y/n): ','s');
if next == 'y',
hold off;
clg
contour(Z,20)
X0=input(' Input guess for Minimum: [x,y]: ');
disp('please wait...');
format bank;
Xmin=fmins('snowmn',X0);
result=[Xmin,snowsl(Xmin(1),Xmin(2))];
disp(' Xmin Ymin Zmin');
disp(result);
X0=input(' Input guess for Maximum: [x,y]: ');
disp('please wait...');
Xmax=fmins('snowmx',X0);
result=[Xmax,snowsl(Xmax(1),Xmax(2))];
disp(' Xmax Ymax Zmax');
disp(result);
end

Draws the Taylor approximation to a sine curve of degree 2k-1

% Draws the Taylor approximation to a sine curve of degree 2k-1
k=input('Type the number of terms of the Taylor series for sine: ');
x=0:.05:2*pi;
z=sin(x);
plot(x,z,'g') %Actual sine curve is drawn in green
hold on
pause % This holds up execution until Enter is pressed
w=x;
y=x;
s=-1;
for j=1:k-1
w=w.*x.*x/(2*j*(2*j+1));
y=y+s*w;
s=-s;
end
plot(x,y,'r') % approximation is drawn in red
hold off

Draws the Taylor approximation to a sine curve or order n

% Draws the Taylor approximation to a sine curve or order n but here the scaling is left fixed after the actual sine curve is drawn

k=input('Type the number of terms of the Taylor series for sine: ');

x=0:.05:6.30;

z=sin(x);
plot(x,z,'g') %Actual sine curve is drawn in green

axis(axis) % Freeze the scaling
hold on

pause

y=x;
s=-1;
w=x;
for j=1:k-1
w=w.*x.*x/(2*j*(2*j+1));
y=y+s*w;
s=-s;
end;
plot(x,y,'r') %Approximation is drawn in red
hold off

vderpol solves the second order Van der Pol equation

function vderpol
%
% vderpol solves the second order Van der Pol equation
%
% w" -k*(1-w*w)*w' + w = 0
%
% as an autonomous first order system in v (=w') and w
% over a time interval 'dt'. It prompts for starting
% coordinates 'w0' and 'v0,' and requests
% permission to complete another time step.
% The derivatives specifying the equation are given in
% vdplfn.m
%
% vderpol.m uses: vdplfn.m, ode23k.m
%

% Initialization
disp('First order system for Van der Pol equation.')
disp(' ')
disp('When you press RETURN you will see an empty plot of v vs w.')
disp('You will be asked for the model parameter k, and initial coordinates.')
disp('Press RETURN to end any pauses.')
pause
more = 'y';
dt = 0.5; trace = 0; tol = 0.05;
xp = [-5 5]; yp = [-5 5];
k = 0;
clg; hold off
plot(xp,yp,'.')
title('1st order system for Van der Pol equation')
xlabel('w')
ylabel('v')
hold on
pause

% Input parameters
k = input(' Model parameter k : ');
while more ~= 'n'
init = input('Initial coordinates (from -4 to 4), enter as [w0 v0]: ');
dt=input('Suggest a time step (try 0.5 if you have no better suggestion): ');
disp('please wait...')
y0 = init'; t0 = 0; tf = dt; step = 'y'; mti=1;
while step == 'y'
%
[tout,yout] = ode23k('vdplfn', k, t0, tf, y0, tol, trace);
%
xp = yout(:,1); yp=yout(:,2);
nti = length(tout); mtf=mti+nti-1;
xpp(mti:mtf)=tout(1:nti);
ypp(mti:mtf)=xp(1:nti); mti=mtf+1;
plot(xp,yp,'r');
pause
disp(' t w(t) v(t)')
disp([tf xp(nti) yp(nti) ])
step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
more = input('Do you want another solution? (y/n): ','s');
end
hold off
clg
disp('Press ENTER for a plot of the last solution: w(t) vs t ');
pause
plot(xpp,ypp,'g');
title('Van der Pol equation')
xlabel('t')
ylabel('w')

Draws a zigzag based on lengths of 100

% Draws a zigzag based on lengths of 100 and l and angles theta1, theta2 (input in degrees).
% This is the most basic version where you have to provide the size of the graphics window and the number of steps.

l=input('Type the value of l ');
theta1=input('Type the angle theta1 in degrees ');
theta2=input('Type the angle theta2 in degrees ');
steps=input('Type the number of steps ');
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; % This is to make the window square

theta1=theta1*pi/180; % Convert to radians!
theta2=theta2*pi/180;



x1=100; y1=0; x2=100+l; y2=0;
% (x2,y2) is the point reached at the end of step 1

x=[0 x1]; y=[0 y1]; plot(x,y,'g');
% This plots the line joining (0,0) to (100,0) in green

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


hold on % So that all plots are made on the same picture

x=[x1 x2]; y=[y1 y2]; plot(x,y,'r');
% This plots the line joining (100,0) to (100+l,0) in red

for j=1:steps-1 % Now we do the zigzag, starting with step 2

x1=x2+100*cos(j*theta1);
y1=y2+100*sin(j*theta1);
% These calculate the end of the next 'zig'

x=[x2 x1]; y=[y2 y1]; plot(x,y,'g');
% This plots the next 'zig' in green

x2=x1+l*cos(j*theta2);
y2=y1+l*sin(j*theta2);
% These calculate the end of the next 'zag'

x=[x1 x2]; y=[y1 y2]; plot(x,y,'r');
% This plots the next 'zag' in red
end;

hold off

Loads a noisy signal similar to fftdemo.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fourdat.
%
% Loads a noisy signal similar to fftdemo.m. Same as that used by
% fouran.m
%
% Usage:-
% fourdat
%
% and press return as required
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo off
% load the data
load noisyt.mat;
load noisys.mat;

disp('Press RETURN for a plot of some of the data')
pause
plot(t(1:50),s(1:50),'g')

load in signal data and then Fourier analyse

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fouran.m
%
% load in signal data and then Fourier analyse
%
%
% Usage:-
% fouran
%
% press return as required
%
% 'fourdata' loads the same data with no Fourier analysis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function fouran

load noisyt.mat
load noisys.mat

%
disp('times and signal values read in OK. No. of sample points is')
nt=length(t)

%make a plot of some of these
np=min(256,nt);nt2=nt/2;np2=np/2;
disp('Press RETURN for a plot of first so many values')
pause
xp=t(1:np); yp=s(1:np);
plot(xp, yp,'g+')
pause

% work out time scale and hence frequency scale
dt=t(2)-t(1); % the time spacing
f=(1:np2-1)/(dt*(np-1)); % the corrresponding frequency values
T=1 ./f; % a corresponding period scale for plotting

% Take the Fourier transform
yf=fft(yp,np); % take the Fourier transform
%
spect=yf(2:np2).*conj(yf(2:np2)); % this is the mod^2 of the
clg % frequency strength
disp(' Now plot against frequency....press RETURN')
pause

% put on a log/linear plot
semilogy(f,spect,'g+')
title(' Frequency spectrum')
xlabel('f /secs')
ylabel('P')
pause

% a log-log plot this time showing the periodicity
disp(' Now plot against period....press RETURN')
pause
loglog(T,spect,'r+')
title(' Periodicity spectrum')
xlabel('T secs')
ylabel('P')

fodesol obtains an approximate solution to the ODE

function fodesol(F,tmin,tmax,xmin,xmax)
%
%
% function fodesol(F,tmin,tmax,xmin,xmax)
%
% fodesol obtains an approximate solution to the ODE
% dx/dt = f(x,t) by the Runge-Kutta-Fehlberg
% method ODE23 supplied with MATLAB (2nd/3rd
% order version) and compares the solutions
% with a "grain" (or field) plot.
%
% F is a string containing the name of a real
% function f(x,t), supplied in an M-file.
% e.g. fnxt.m
% tmin,tmax and
% xmin,xmax are limits for the plot
%
% usage e.g.
% fodesol('fnxt',0, 2.5, -3 , 3)
%
% fodesol.m uses: odegr.m and mode23.m
%

% Initialization
disp('Here is a grain or field plot of dx/dt = f(x,t)')
disp('Press Enter to see it')
pause
more = 'y';
clg; hold off

% Grain plot
%axis([tmin tmax xmin xmax])
odegr(F,tmin,tmax,xmin,xmax)
hold on
axis(axis)
% Superimpose numerical solution through given point
while more == 'y'
t0 = input('Pick initial t-value: ');
x0 = input('Pick initial x-value: ');
disp('please wait...')
[xp,yp] = mode23(F,t0,tmax,x0);
% use a spline to make it look smoother
n=length(xp); x1=xp(1); x2=xp(n);
dx=(x2-x1)/50; xi=x1:dx:x2;
yi=spline(xp',yp',xi);
plot(xi,yi,'r'); pause
more = input('Do you want another solution? (y/n): ','s');
end
disp(' Initial and final points (t,x) of solution were')
l=length(yp);
tf=xp(l); xf=yp(l);
disp(' Initial Final')
disp([t0 tf])
disp([x0,xf])
hold off

Flu epidemic data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fludat.m
%
% Flu epidemic data
%
% Usage:
% fludat
%
% Loads in matrix alldat(ndata,2)
% and vectors xi (tdays), yi (I)
% which are its columns
% tdays is the time in days
% I is the number of infectives (ill boys) at time t
% The number of data points (ndata) is also returned
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo off;
hold off;

%Data goes in here
% tdays I
alldat=[
3 24
4 79
5 230
6 301
7 256
8 245
9 201
10 138
11 78
12 25
13 12
14 9];

% end of the data
% define a few convenient variables
tdays = alldat(:,1); xi=tdays;
I = alldat(:,2); yi=I;
disp('No. of data points:')
ndata = length(tdays);
disp(ndata)

%plot the data
disp(' Press RETURN for plot')
pause
%define limits of the plot (bottom LH, top RH)
xp=[0 16]; yp=[0 800];
plot(xp,yp,'.w')
hold on

plot(xi,yi,'*g')
title('Flu infections')
xlabel('days')
ylabel('I')
hold off

Estimates the convergence rate of a given sequence

function Rate_k = name(Root, Id);

% f_rate.m

%

% Estimates the convergence rate of a given sequence

% of scalars (Root=Vector) or vectors (Root=Matrix).

%

% Usage: K = f_rate ( Root, id)

% ==or== K = f_rate ( Root )

%

% where Root is of size (m x n), with n = dimension, and

% id=1 is to show all intermediate estimates.

% For usage 2, id=0;

% in all cases the very last estimate is shown.

% For example, to check the convergence of the iterates,

% (1, 2, 3, 4, 4.5, 5, 5.11)

% USE Root = [1 2 3 4 4.5 5 5.11]';

% kc_1 = f_rate(Root,1) %% Giving 1.3914 ... 2.8614

% kc_2 = f_rate(Root) %% Giving 1.7101 and 2.8614



if (nargin < 1), help f_rate; return, end

if nargin==1, Id=0; end

[m n] = size(Root) ; % Check if m is more than 4

if (m <>= 4), Root=Root'; [m n]=size(Root); end

abort_id = 0;

x_exact = Root(m,:) ; % Take as exact solution "x"

keep = [] ; Rate_k = []; i=3;

if m >= 4

while (i <> 3 or m = 3

d_n_m = norm(Root(i-2,:)-x_exact) ; % d_(n-1)

d_n = norm(Root(i-1,:)-x_exact) ; % d_n

d_n_p = norm(Root(i,:) -x_exact) ; % d_(n+1)

if (d_n_m < eps) d_n_m = eps*7; end % In case of 0

if (d_n < eps) d_n = eps*7; end % In case of 0

if (d_n_m < eps) d_n_m = eps*7; end % In case of 0

if (d_n < eps) d_n = eps*7; end % In case of 0

if (d_n_p < eps)

d_n_p = eps*7 ; % In case of 0

i = m - 1; abort_id = 1; % No need to continue

end

top = d_n_p / d_n ; % numerator

bottom = d_n / d_n_m ; % denominator

if (top == 1) top = 1 + eps*7 ; end % Safeguard for Log

if (bottom == 1) bottom = 1 + eps*7 ; end

rate_k = log(top) / log(bottom) ; % Estimate of the order

if (abort_id == 0 & rate_k > 0), keep = [keep; rate_k]; end

i = i + 1;

end % End WHILE-loop

disp(' ');

if Id == 1

Rate_k = keep ;

disp(' (Individual estimates for rate)') ;

else

m = length(keep); m=keep(m);

fprintf('Average Rate =%5.1f\n',mean(keep));

fprintf('Last estimate=%5.1f\n',m);

Rate_k = [mean(keep) m];

end % End for IF

else

fprintf('Rate_k is not estimated, because size = %d\n',m)

disp('Not enough iterates have been supplied for f_rate.m');

help f_rate ;

end; % END the main IF m>4

M-file to plot the Fibonacci numbers f(i) less than 1000

%M-file to plot the Fibonacci numbers f(i) less than 1000.
%The Fibonacci numbers are defined by f(1)=f(2)=1
%and f(i+2)=f(i+1)+f(i) for i=1,2,3,...
%They are calculated and then plotted with i along the
%horizontal axis and f(i) along the vertical axis.

f=[1 1];
k=1;
while f(k) < 1000
f(k+2)=f(k+1)+f(k);
k=k+1;
end
f
plot(f)

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)