Chapter 9 M-Files

M-File Home

Script Files

Function Files

ShowTraj

% Script File: ShowTraj
% Plots a family of solutions to y'(t) = -5y(t)

close all
t= linspace(-.1,.4);
y = exp(-5*t);
plot(t,y);
axis([-.1 .4 0 2])
hold on
for c=linspace(0,4,21)
   plot(t,c*y)
end
plot(0,1,'o')
plot(t,y,'LineWidth',2)
hold off
title('Solutions to y''(t)  =  -5 y(t)')

ShowEuler

% Script File: ShowEuler
% Plots a family of solutions to y'(t) = -5y(t)

close all
t= linspace(-.1,.4);
y = exp(-5*t);
plot(t,y);
axis([-.1 .4 0 2])
hold on
tc = 0;
yc = 1;
plot(tc,yc,'*')

hold on
for k=0:4
   title(sprintf('k=%1.0f,  Click t(%1.0f)',k,k+1))
   [tnew,z] = ginput(1);
   hc = tnew-tc;
   fc = -5*yc;
   ynew = yc + hc*fc;
   plot([tc tnew],[yc ynew],'--',tnew,ynew,'o')
   tc = tnew;
   yc = ynew;  
end
title('Five Steps of Euler Method (y''=-5y, y(0)=1) ')
hold off

ShowFixedEuler

% Script File: ShowFixedEuler
% Plots the error in a fixed h euler solution to
% y' = -y, y(0) = 1,  on [0,5]

close all
tol=.01;
[tvals,yvals] = FixedEuler('f1',1,0,5,1,tol);
plot(tvals,exp(-tvals)-yvals)
title('Fixed h Euler Error for y''=-y, y(0) = 1')
xlabel(sprintf('tol = %5.3f,   n = %4.0f',tol,length(tvals)))

ShowTrunc

% Script File: ShowTrunc
% Illustrates local truncation error

close all
a=-5;
h=.1;
y0=1;
t = linspace(-h,4*h,100);
clf
y = exp(a*t);
plot(t,y)
text(.2,1.18,sprintf('y'' = %4.1fy, y(0)  = 1',a))
text(.2,1.1,'* = exact solution')
text(.2,1.02,'o = computed solution')
hold on
plot(0,y0,'*')
y1 = (1 + a*h)*y0;
y = (y1/exp(a*h))*exp(a*t);
plot(t,y)
plot(h,exp(a*h),'*')
plot(h,(y1/exp(a*h))*exp(a*h),'o')
pause(1)
y2 = (1 + a*h)*y1;
y = (y2/exp(2*a*h))*exp(a*t);
plot(t,y)
plot(2*h,exp(a*2*h),'*')
plot(2*h,(y2/exp(2*a*h))*exp(a*2*h),'o')
y3 = (1 + a*h)*y2;
y = (y3/exp(3*a*h))*exp(a*t);
plot(t,y)
plot(3*h,exp(a*3*h),'*')
plot(3*h,(y3/exp(3*a*h))*exp(a*3*h),'o')
hold off

%
title('Euler Solution (h=0.1) of y''=-5y, y(0) = 1')

EulerRoundoff

% Script File: EulerRoundOff
% Global error for the fixed step Euler method in 3-digit floating point arithmetic
% applied to y' = -y across [0,1].

close all
t = linspace(0,1);
e = exp(-t);
for n = 140:20:180
   tvals = linspace(0,1,n+1);
   yvals = 0*tvals;
   one = represent(1);
   yc = one;
   factor = float(one,float(one,represent(n),'/'),'-');
   yvals(1) = 1;
   for k=1:n
      yc = float(factor,yc,'*');
      yvals(k+1) = convert(yc);
   end
   plot(tvals,abs(exp(-tvals)-yvals))
   hold on
   axis([0 1 0 .05])	
end
hold off
title('Euler in 3-Digit Arithmetic')
text(.76,.005,'h = 1/140')
text(.62,.012,'h = 1/160')
text(.42,.02,'h = 1/180')

ShowAB

% Script File: ShowAB
% Plots absolute error for fixed step size Adams-Bashforth
% solution to y' = y, y(0) = 1 across [0,5].

close all
E = zeros(5,5);
for i=1:5
   n = 16*2^i;
   Exact = exp(-linspace(0,5,n+1)');
   for k=1:5  
      [tvals,yvals] = FixedAB('f1',0,1,5/n,k,n);
      E(i,k) = max(abs(yvals-Exact));
   end
end
semilogy([32 64 128 256 512]',E)
title('Adams-Bashforth on y''(t) = -y(t), y(0) = 1, 0<=t<=5, h = 5/n')
xlabel('n (Number of Steps)')
ylabel('Maximum absolute error')
text(530,E(5,1),'k=1')
text(530,E(5,2),'k=2')
text(530,E(5,3),'k=3')
text(530,E(5,4),'k=4')
text(530,E(5,5),'k=5')
axis([0 600 1e-12 1e-1])

ShowPC

% Script File: ShowPC
% Plots absolute error for fixed step size Adams PC
% solution to y' = y, y(0) = 1 across [0,5].

close all
E = zeros(5,5);
for i=1:5
   n = 16*2^i;
   Exact = exp(-linspace(0,5,n+1)');
   for k=1:5  
      [tvals,yvals] = FixedPC('f1',0,1,5/n,k,n);
      E(i,k) = max(abs(yvals-Exact));
   end
end
semilogy([32 64 128 256 512]',E)
title('Adams PC on y''(t) = -y(t), y(0) = 1, 0<=t<=5, h = 5/n')
xlabel('n (Number of Steps)')
ylabel('Maximum absolute error')
text(530,E(5,1),'k=1')
text(530,E(5,2),'k=2')
text(530,E(5,3),'k=3')
text(530,E(5,4),'k=4')
text(530,E(5,5),'k=5')
axis([0 600 1e-14 1e-1])

ShowRK

% Script File: ShowRK
% Plots absolute error for fixed step size Runge-Kutta
% solution to y' = y, y(0) = 1 across [0,5].

close all
E = zeros(5,5);
for i=1:5
   n = 16*2^i;
   Exact = exp(-linspace(0,5,n+1)');
   for k=1:5  
      [tvals,yvals] = FixedRK('f1',0,1,5/n,k,n);
      E(i,k) = max(abs(yvals-Exact));
   end
end
semilogy([32 64 128 256 512]',E)
title('Runge-Kutta on y''(t) = -y(t), y(0) = 1, 0<=t<=5, h = 5/n')
xlabel('n (Number of Steps)')
ylabel('Maximum absolute error')
text(530,E(5,1)+.003,'k=1')
text(530,E(5,2),'k=2')
text(530,E(5,3),'k=3')
text(530,E(5,4),'k=4')
text(530,E(5,5),'k=5')

ShowKepler

% Script File: ShowKepler
% Applies ODE23 and ODE45 to a system of differential equations
% that define an elliptical orbit.

close all
clc

% A simple call to ode23.

tInitial = 0;
tFinal   =  2*pi;
uInitial = [ .4; 0 ; 0 ; 2];
tSpan = [tInitial tFinal];
[t, u] = ode23('Kepler', tSpan, uInitial);
nSteps = length(t)-1;
plot(u(:,1),u(:,3))
axis('equal')
title('Kepler Problem: ode23 with Default Tolerances')
xlabel(sprintf('Number of Steps = %5d',nSteps))
figure
plot(t(2:length(t)),diff(t))
title('Kepler Problem: ode23 with Default Tolerances')
ylabel('Step Length')
xlabel('t')
figure
subplot(2,2,1), plot(t,u(:,1)), title('x(t)')
subplot(2,2,2), plot(t,u(:,3)), title('y(t)')
subplot(2,2,3), plot(t,u(:,3)), title('x''(t)')
subplot(2,2,4), plot(t,u(:,4)), title('y''(t)')

% A call with specified output times.

tSpan = linspace(tInitial,tFinal,20);
[t, u] = ode23('Kepler', tSpan, uInitial);
xvals = spline(t,u(:,1),linspace(0,2*pi));
yvals = spline(t,u(:,3),linspace(0,2*pi));
figure
plot(xvals,yvals,u(:,1),u(:,3),'o')
axis('equal')
title('Kepler Problem: ode23 with Specified Output Times')
legend('Spline Fit','ode23 Output',0)

% A call with a more stringent tolerances

tSpan = [tInitial tFinal];
options = odeset('AbsTol',.00000001,'RelTol',.000001,'stats','on');
disp(sprintf('\n Stats for ODE23 Call:\n'))
[t, u] = ode23('Kepler', tSpan, uInitial,options);
nSteps = length(t)-1;
figure
plot(u(:,1),u(:,3))
axis('equal')
title('Kepler Problem: ode23 with RelTol = 10^{-6} and AbsTol = 10^{-8}')
xlabel(sprintf('Number of Steps = %5d',nSteps))
figure
plot(t(2:length(t)),diff(t))
title('Kepler Problem: ode23 with RelTol = 10^{-6} and AbsTol = 10^{-8}')
ylabel('Step Length')
xlabel('t')


% Use ODE45 on the same problem.

tSpan = [tInitial tFinal];
options = odeset('AbsTol',.00000001,'RelTol',.000001,'stats','on');
disp(sprintf('\n Stats for ode45 Call:\n'))
[t, u] = ode45('Kepler', tSpan, uInitial,options);
nSteps = length(t)-1;
figure
plot(u(:,1),u(:,3))
axis('equal')
title('Kepler Problem: ode45 with RelTol = 10^{-6} and AbsTol = 10^{-8}')
xlabel(sprintf('Number of Steps = %5d',nSteps))
figure
plot(t(2:length(t)),diff(t))
title('Kepler Problem: ode45 with RelTol = 10^{-6} and AbsTol = 10^{-8}')
ylabel('Step Length')
xlabel('t')

FixedEuler

  function [tvals,yvals] = FixedEuler(fname,y0,t0,tmax,M2,tol)
% [tvals,yvals] = FixedEuler(fname,y0,t0,tmax,M2,tol)
% Fixed step Euler method.
%
% fname is a string that names a function of the form f(t,y). 
% M2 a bound on the second derivative of the solution to
%                 y' = f(t,y),   y(t0) = y0
% on the interval [t0,tmax].
%
% Determine positive n so that if tvals = linspace(t0,tmax,n), then
% y(i) is within tol of the true solution y(tvals(i)) for i=1:n.
 
n = ceil(((tmax-t0)^2*M2)/(2*tol))+1;
h = (tmax-t0)/(n-1);
yvals = zeros(n,1);
tvals = linspace(t0,tmax,n)';
yvals(1) = y0;
for n=1:n-1
   fval = feval(fname,tvals(n),yvals(n));
   yvals(n+1) = yvals(n)+h*fval;
end

ABStart

  function [tvals,yvals,fvals] = ABStart(fname,t0,y0,h,k)
% [tvals,yvals,fvals] = ABStart(fname,t0,y0,h,k)
% Generates enough starting values for the kth order Adams-Bashforth method.
% Uses k-th order Runge-Kutta to generate approximate solutions to
%
%               y'(t) = f(t,y(t))   y(t0) = y0
%
% at t = t0, t0+h, ... , t0 + (k-1)h. 
%
% fname is a string that names the function f.
% t0 is the initial time.
% y0 is the initial value.
% h is the step size.
% k is the order of the RK method used.
%
% tvals is a column vector with tvals(j) = t0 + (j-1)h, j=1:k
% yvals is a matrix with yvals(j,:) = approximate solution at tvals(j), j=1:k
% For j =1:k, fvals(:,j) = f(tvals(j),yvals(j,:)).

tc = t0;    
yc = y0;     
fc = feval(fname,tc,yc);
tvals = tc; 
yvals = yc'; 
fvals = fc;
for j=1:k-1
   [tc,yc,fc] = RKstep(fname,tc,yc,fc,h,k);
   tvals = [tvals; tc]; 
   yvals = [yvals; yc']; 
   fvals = [fc fvals];
end

ABStep

  function [tnew,ynew,fnew] = ABStep(fname,tc,yc,fvals,h,k) 
% [tnew,ynew,fnew] = ABStep(fname,tc,yc,fvals,h,k)  
% Single step of the kth order Adams-Bashforth method.
%
% fname is a string that names a function of the form f(t,y)
% where t is a scalar and y is a column d-vector.
%
% yc is an approximate solution to y'(t) = f(t,y(t)) at t=tc.
%
% fvals is an d-by-k matrix where fvals(:,i) is an approximation
% to f(t,y) at t = tc +(1-i)h, i=1:k
%
% h is the time step. 
%
% k is the order of the AB method used, 1<=k<=5.
%
% tnew=tc+h.
% ynew is an approximate solution at t=tnew.
% fnew = f(tnew,ynew).

if k==1
   ynew = yc + h*fvals;
elseif k==2
   ynew = yc + (h/2)*(fvals*[3;-1]);
elseif k==3
   ynew = yc + (h/12)*(fvals*[23;-16;5]);
elseif k==4
   ynew = yc + (h/24)*(fvals*[55;-59;37;-9]);
elseif k==5
   ynew = yc + (h/720)*(fvals*[1901;-2774;2616;-1274;251]);
end
tnew = tc+h;
fnew = feval(fname,tnew,ynew);

FixedAB

  function [tvals,yvals] = FixedAB(fname,t0,y0,h,k,n)
% [tvals,yvals] = FixedAB(fname,t0,y0,h,k,n)
% Produces an approximate solution to the initial value problem
%
%      y'(t) = f(t,y(t))     y(t0) = y0
%
% using a strategy that is based upon a k-th order
% Adams-Bashforth method. Stepsize is fixed.
%
% fname = string that names the function f.
% t0 = initial time.
% y0 = initial condition vector.
% h  = stepsize.
% k  = order of method. (1<=k<=5).
% n  = number of steps to be taken,
%
% tvals is a column vector with tvals(j) = t0 + (j-1)h, j=1:n+1
% yvals is a matrix with yvals(j,:) = approximate solution at tvals(j), j=1:n+1

[tvals,yvals,fvals] = ABStart(fname,t0,y0,h,k);
tc = tvals(k);
yc = yvals(k,:)';
fc = fvals(:,k);

for j=k:n
   % Take a step and then update.
   [tc,yc,fc] = ABstep(fname,tc,yc,fvals,h,k);
   tvals = [tvals; tc]; 
   yvals = [yvals; yc'];
   fvals = [fc fvals(:,1:k-1)];
end

AMStep

  function [tnew,ynew,fnew] = AMstep(fname,tc,yc,fvals,h,k) 
% [tnew,ynew,fnew] = AMstep(fname,tc,yc,fvals,h,k)
% Single step of the kth order Adams-Moulton method.
%
% fname is a string that names a function of the form f(t,y)
% where t is a scalar and y is a column d-vector.
%
% yc is an approximate solution to y'(t) = f(t,y(t)) at t=tc.
%
% fvals is an d-by-k matrix where fvals(:,i) is an approximation
% to f(t,y) at t = tc +(2-i)h, i=1:k
%
% h is the time step. 
%
% k is the order of the AM method used, 1<=k<=5.
%
% tnew=tc+h
% ynew is an approximate solution at t=tnew
% fnew = f(tnew,ynew).


if k==1
   ynew = yc + h*fvals;
elseif k==2
   ynew = yc + (h/2)*(fvals*[1;1]);
elseif k==3
   ynew = yc + (h/12)*(fvals*[5;8;-1]);
elseif k==4
   ynew = yc + (h/24)*(fvals*[9;19;-5;1]);
elseif k==5
   ynew = yc + (h/720)*(fvals*[251;646;-264;106;-19]);
end
tnew = tc+h;
fnew = feval(fname,tnew,ynew);

PCStep

  function [tnew,yPred,fPred,yCorr,fCorr] = PCStep(fname,tc,yc,fvals,h,k)  
% [tnew,yPred,fPred,yCorr,fCorr] = PCStep(fname,tc,yc,fvals,h,k)
% Single step using the kth-order Adams predictor-corrector framework.
%
% fname is a string that names a function of the form f(t,y)
% where t is a scalar and y is a column d-vector.
%
% yc is an approximate solution to y'(t) = f(t,y(t)) at t=tc.
%
% fvals is an d-by-k matrix where fvals(:,i) is an approximation
% to f(t,y) at t = tc +(1-i)h, i=1:k
%
% h is the time step. 
%
% k is the order of the Adams methods used, 1<=k<=5.
%
% tnew=tc+h, 
% yPred is the predicted solution at t=tnew
% fPred = f(tnew,yPred)
% yCorr is the corrected solution at t=tnew
% fCorr = f(tnew,yCorr).

[tnew,yPred,fPred] = ABstep(fname,tc,yc,fvals,h,k);  
[tnew,yCorr,fCorr] = AMstep(fname,tc,yc,[fPred fvals(:,1:k-1)],h,k);

FixedPC

  function [tvals,yvals] = FixedPC(fname,t0,y0,h,k,n)
% [tvals,yvals] = FixedPC(fname,t0,y0,h,k,n)
% Produces an approximate solution to the initial value problem
%
%      y'(t) = f(t,y(t))     y(t0) = y0
%
% using a strategy that is based upon a k-th order
% Adams Predictor-Corrector framework. Stepsize is fixed.
%
% fname = string that names the function f.
% t0 = initial time.
% y0 = initial condition vector.
% h  = stepsize.
% k  = order of method. (1<=k<=5).
% n  = number of steps to be taken,
%
% tvals is a column vector with tvals(j) = t0 + (j-1)h, j=1:n+1
% yvals is a matrix with yvals(j,:) = approximate solution at tvals(j), j=1:n+1.

[tvals,yvals,fvals] = ABStart(fname,t0,y0,h,k);
tc = tvals(k);
yc = yvals(k,:)';
fc = fvals(:,k);

for j=k:n
   % Take a step and then update.
   [tc,yPred,fPred,yc,fc] = PCstep(fname,tc,yc,fvals,h,k);
   tvals = [tvals; tc];
   yvals = [yvals; yc'];
   fvals = [fc fvals(:,1:k-1)];
end

RKStep

  function [tnew,ynew,fnew] = RKStep(fname,tc,yc,fc,h,k)
% [tnew,ynew,fnew] = RKStep(fname,tc,yc,fc,h,k)
% Single step of  the kth order Runge-Kutta method.
%
% fname is a string that names a function of the form f(t,y)
% where t is a scalar and y is a column d-vector.
%
% yc is an approximate solution to y'(t) = f(t,y(t)) at t=tc.
%
% fc = f(tc,yc).
%
% h is the time step. 
%
% k is the order of the Runge-Kutta method used, 1<=k<=5.
%
% tnew=tc+h, ynew is an approximate solution at t=tnew, and 
% fnew = f(tnew,ynew).

if k==1
   k1 = h*fc;
   ynew = yc + k1;
elseif k==2
   k1 = h*fc;
   k2 = h*feval(fname,tc+h,yc+k1);
   ynew  = yc + (k1 + k2)/2;
elseif k==3 
   k1 = h*fc;
   k2 = h*feval(fname,tc+(h/2),yc+(k1/2));
   k3 = h*feval(fname,tc+h,yc-k1+2*k2);
   ynew  = yc + (k1 + 4*k2 + k3)/6;
elseif k==4
   k1 = h*fc;
   k2 = h*feval(fname,tc+(h/2),yc+(k1/2));   
   k3 = h*feval(fname,tc+(h/2),yc+(k2/2));  
   k4 = h*feval(fname,tc+h,yc+k3);  
   ynew  = yc + (k1 + 2*k2 + 2*k3 + k4)/6;
elseif k==5
   k1 = h*fc;
   k2 = h*feval(fname,tc+(h/4),yc+(k1/4));   
   k3 = h*feval(fname,tc+(3*h/8),yc+(3/32)*k1+(9/32)*k2);    
   k4 = h*feval(fname,tc+(12/13)*h,yc+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3);               
   k5 = h*feval(fname,tc+h,yc+(439/216)*k1 - 8*k2 + (3680/513)*k3 -(845/4104)*k4);
   k6 = h*feval(fname,tc+(1/2)*h,yc-(8/27)*k1 + 2*k2 -(3544/2565)*k3 + (1859/4104)*k4 - (11/40)*k5);
   ynew  = yc + (16/135)*k1 + (6656/12825)*k3 + (28561/56430)*k4 - (9/50)*k5 + (2/55)*k6;
end
tnew = tc+h;
fnew = feval(fname,tnew,ynew);

FixedRK

  function [tvals,yvals] = FixedRK(fname,t0,y0,h,k,n)
% [tvals,yvals] = FixedRK(fname,t0,y0,h,k,n)
% Produces approximate solution to the initial value problem
%
%      y'(t) = f(t,y(t))     y(t0) = y0
%
% using a strategy that is based upon a k-th order
% Runge-Kutta method. Stepsize is fixed.
%
% fname = string that names the function f.
% t0 = initial time.
% y0 = initial condition vector.
% h  = stepsize.
% k  = order of method. (1<=k<=5).
% n  = number of steps to be taken,
%
% tvals is a column vector with tvals(j) = t0 + (j-1)h, j=1:n+1
% yvals is a matrix with yvals(j,:) = approximate solution at tvals(j), j=1:n+1

tc = t0;
yc = y0;
tvals = tc;
yvals = yc';
fc = feval(fname,tc,yc);
for j=1:n
   [tc,yc,fc] = RKstep(fname,tc,yc,fc,h,k);
   yvals = [yvals; yc' ];
   tvals = [tvals; tc];
end

Kepler

  function up = Kepler(t,u)
% up = Kepler(t,u)
% t (time) is a scalar and u is a 4-vector whose components satisfy
%
%            u(1) = x(t)      u(2) = (d/dt)x(t)
%            u(3) = y(t)      u(4) = (d/dt)y(t)
%
% where (x(t),y(t)) are the equations of motion in the 2-body problem.
%
% up is a 4-vector that is the derivative of u at time t.

r3 = (u(1)^2 + u(3)^2)^1.5;
up = [  u(2)     ;...
       -u(1)/r3  ;...
        u(4)     ;...
       -u(3)/r3] ;

f1

function yp = f1(t,y)
yp = -y;