Chapter 9: Problem Solutions

Problems Home

P9.1.1 P9.2.1 P9.3.1
P9.1.2 P9.2.2 P9.3.2
P9.1.3 P9.2.3 P9.3.3
P9.1.4 P9.2.4
P9.2.5
P9.2.6
P9.2.7

 


P9.1.1

% Problem P9_1_1
%
% Examine Euler's method

close all
[tvals yvals] = FixedEuler('MyF911',1,1,4,10,.01');
plot(tvals,yvals)

  function yp = MyF911(t,y)
yp = -t*y + 1/y^2;

P9.1.2

% Problem P9_1_2
% 
% Apply Euler's method to a system

clc
[T, Y] = ODE45('MyF912',linspace(0,1,20), [3;-1;4]);
[m,n] = size(Y);
yat1 = Y(m,n);
error = 1;
n = 10;

disp('    n     y(1)     error')
disp('--------------------------')
while error >=.001
   n = 2*n;
   t = linspace(0,1,n+1)';
   h = 1/n;
   u = [3;-1;4];
   for k=1:n
      u = u + h*MyF912(t(k),u);
   end
   error = abs(yat1 - u(3))/abs(yat1);
   disp(sprintf(' %4.0f  %7.4f   %7.4f',n,u(3),error))
end

  function up = MyF912(t,u)
x  = u(1);
xp = u(2);
y  = u(3);
up = zeros(3,1);
up(1) = xp;
up(2) = (3 - sin(t))*xp + x/(1 + y^2);
up(3) = -cos(t)*y - xp/(1+t^2);

P9.1.3

% Problem  P9_1_3
% 
% Apply Euler's method to a system.

n = 100;
t = linspace(0,3,n+1)';
h = 3/n;
Y = zeros(2,n+1);
Y(:,1) = [2;-1];
for k=1:n
   Y(:,k+1) = Y(:,k) + h*MyF913(t(k),Y(:,k));
end
plot(t,Y(1,:),t,Y(2,:))


  function yp = MyF913(t,y)
yp = [-1 4; -4 -1]*y;

P9.1.4

Not Available

P9.2.1

% Problem 9_2_1
%
% Examine RKF45 applied to y' = t^d, y(1) = 1 which has solution 
% y(t) = (t^(d+1) + d)/(d+1). Note that the d+2nd derivative of the solution
% is zero suggesting that a d+1st order method will be exact.

tinitial = 1;
tfinal = 2;
yinitial = 1;
h = tfinal-tinitial;
global d
clc 
disp(' d      ynew error       znew error')
disp('-------------------------------------')
for d=0:5
   yexact = (tfinal^(d+1) + d)/(d+1);
   [tnew,ynew,znew] = RKF45step('MyF921',tinitial,yinitial,h);
   disp(sprintf('%1.0f      %8.6e     %8.6e',d,abs(yexact-ynew), abs(yexact-znew)))
end

  function yp = MyF921(t,y)
global d
yp = t^d;

P9.2.2

% Problem P9_2_2
%
% Error at t= 1 for fixed stepsize RK applied to y' = -y , y(0) = 1.

clc
fevals = zeros(1,5);
guess = [1000 250 26 8 4];
for k=1:5
   n=guess(k);
   while n>0 
      h = 1/n;
      tc = 0; yc = 1; fc = -1;
      for j=1:n
         [tc,yc,fc] = RKStep('f1',tc,yc,fc,h,k);
      end
      error = abs(exp(-1)-yc);
      fevals = k*n;
      if k==5
         fevals = 6*k;
      end
      disp(sprintf('k = %1.0f, n = %4.0f, Error = %12.8f, f-evals = %4.0f',k,n,error,fevals))		
      n = input('Enter zero to quit or another positive choice for n:');   
   end	  
end



P9.2.3

% Script 9_2_3
%
% Illustrate ode23 and ode45 on the Apollo orbit problem
global fcount
tinitial = 0;
uinitial = [1.2;0;0;-1.0493575];
tfinal = 6.19216933;
mu = 1/82.45;
close all
clc home
disp(' tol        ode23 steps    f-evaluations')
disp('----------------------------------------')
for tol = [.01 .001   .0001]
   fcount = 0;
   [t,u] = Ode23('MyF923',tinitial,tfinal,uinitial,tol);
   figure
   plot(u(:,1),u(:,3),1-mu,0,'o',-mu,0,'o')
   title(sprintf('Ode23 (tol = %6.1e)',tol))
   disp(sprintf('%6.1e     %5.0f           %5.0f',tol,length(t),fcount))
end

disp(' ')
disp(' ')
disp(' tol        ode45 steps    f-evaluations')
disp('----------------------------------------')
for tol = [.01 .001   .0001  .00001]
   fcount = 0;
   [t,u] = Ode45('MyF923',tinitial,tfinal,uinitial,tol);
   figure
   plot(u(:,1),u(:,3),1-mu,0,'o',-mu,0,'o')
   title(sprintf('Ode45 (tol = %6.1e)',tol))
   disp(sprintf('%6.1e     %5.0f           %5.0f',tol,length(t),fcount))
end


  function up = MyF923(t,u)
%
% t is a real scalar
% u is a column 4-vector
%
% up column 4-vector that reflects the rate of change of u at t. 

global fcount
fcount = fcount + 1;
x  = u(1);
xp = u(2);
y  = u(3);
yp = u(4);
mu = 1/82.5;
mu_star = 1-mu;
r1 = sqrt((x+mu)^2      + y^2);
r2 = sqrt((x-mu_star)^2 + y^2);
up = zeros(4,1);
up(1) = u(2);
up(2) =  2*yp + x - mu_star*(x+mu)/r1^3  - mu*(x-mu_star)/r2^3;
up(3) = u(4);
up(4) = -2*xp + y - mu_star*y/r1^3       - mu*y/r2^3;

P9.2.4

% Problem 9_2_4
% Solve the problem  
%                    y'' = -2y' - 2y + 6          y(0) = 3, y(5) = 3+exp(-5)sin(5)
%
% This has solution y(t) = 3 + exp(-t)sin(t).

% Let u1(t) = y(t) and u2(t) = y'(t) then the problem reformulates as
%
%        u1'(t) = u2(t)                       u1(0) = 3
%        u2'(t) = -2u2(t) - 2u1(t) + 6        u1(5) = 3+exp(-5)sin(5)
%
% Our goal is to solve for u1(t). Let g(mu) be the estimate of u1(5) obtained
% by applying ode23 to the IVP

%        u1'(t) = u2(t)                       u1(0) = 3
%        u2'(t) = -2u2(t) - 2u1(t) + 6        u2(0) = mu
%
% Use fzero to find a zero of G(mu) = g(mu) - u1(5). (Observe that mu = 1 is a root
% since for the true solution u2(0) = y'(0) = 1.)

% We know that mu = 1 is the correct solution so let's start with an initial
% guess mu = .5.

mu_star = fzero('g924fzero',.5);

% Generate the approximate solution and compare.

tspan = linspace(0,5,21)';
[t,u] = ode23('f924ode',tspan,[3;mu_star]);
tau = linspace(0,5,200)';
plot(tau,3+exp(-tau).*sin(tau),t,u(:,1),'o')
legend('3 + exp(-t)sin(t)','Computed approximation',0)


function ypp = f924(t,y,yp)
% The given "second derivative function"  y" = f(t,y,y').

ypp = -2*yp - 2*y + 6;



  function up = f924ode(t,u)
%
% up(1) = u(2);
% up(2) = f924(t,u(1),u(2))

up = [u(2); f924(t,u(1),u(2))];


  function z = g924(mu)
% z is an estimate of u1(5) where
%
%        u1'(t) = u2(t)                       u1(0) = 3
%        u2'(t) = -2u2(t) - 2u1(t) + 6        u2(0) = mu

[t,zvals] = ode23('f924ode',[0 5],[3;mu]);

z = zvals(length(t),1);


function r = g924fzero(mu)
% This will be zero when the initial condition u2(0) = mu
% renders the right value for u1(t) at t = 5.

  r = g924(mu) - (3+exp(-5)*sin(5));
  
  




P9.2.5

Not Available

P9.2.6

% Problem P9_2_6

randn('state',0)
n = 8;

% Generate an n-by-n positive definite matrix.

X = randn(n,n);
A = X'*X;

% Compute its Cholesky factorization via cholSax from Chapter 6...

G = chol(A);  % A = GG', thus to solve Az = b, execute z = G'\(G\b);

% We'll solve y' = Cy, y(0) = ones(n,1) where C is the matrix

C = randn(n,n);
C = C-C';

m = 101;
T = 10;

% Solve the IVP. Get solutions at linspace(0,T,m).

tspan = linspace(0,T,m)';
[t,y] = ode23('F926a',tspan,ones(n,1),[],C);

% Note that y(i,:)' is an approximation of the true solution at t(i).

% Let f(t) = y'(t)*A^{-1}y(t). Compute f(i) = y(t(i)), i=1:m

f = zeros(m,1);
for i=1:m
   f(i) = y(i,:)*(G\(G'\y(i,:)'));
end

% Generate the spline interpolant and integrate

S = spline(t,f);
Q = quad('F926b',0,T,[],[],S)


function yp = F926a(t,y,flag,C)
yp = C*y;

function s = F926b(t,S)
s = ppval(S,t);

P9.2.7

Noty Available

P9.3.1

% Problem P9_3_1
%
% Specialize the Adams methods for the case y' = Ay.

close all
global A
A = [ -3 1 -2; -1 -5 3; 1 3 -4];


y0 = [1;1;1];
n = 100;
h = 2/n;
t0 = 0;
% Solve y' = [ -3 1 -2; -1 -5 3; 1 3 -4]y , y(0) = [1;1;1]
% across the interval [0,2] using h = 1/50.
% Exact solution usingthe matrix exponential.
F = expm(h*A);
tvals = linspace(t0,t0+n*h,n+1);
yexact = y0;
for k=1:n
   yexact = [yexact F*yexact(:,k)];
end

% The Adams-Bashforth solutions.
for k=1:5  
   [tvals,yvals] = AFixedAB(A,t0,y0,h,k,n);
   figure
   E = abs(yexact - yvals');
   plot(tvals,E(1,:),tvals,E(2,:),tvals,E(3,:))
   title(sprintf('AB%1.0f: Error in the solution components',k))
   legend('computed y_{1}(t) error', 'computed y_{2}(t) error','computed y_{3}(t) error')
end

% The Adams-Moulton solutions.
for k=1:5  
   [tvals,yvals] = AFixedAM(A,t0,y0,h,k,n);
   figure
   E = abs(yexact - yvals');
   plot(tvals,E(1,:),tvals,E(2,:),tvals,E(3,:))
   title(sprintf('AM%1.0f: Error in the solution components',k)) 
   legend('computed y_{1}(t) error', 'computed y_{2}(t) error','computed y_{3}(t) error')
end



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

tvals = linspace(t0,t0+h*n,n+1);
yvals = zeros(n+1,length(y0));
[tvals(1:k),yvals(1:k,:)] = ABStart('MyF931',t0,y0,h,k);

% Store the solutions as column vectors to make the iteration simpler.

yvals = yvals'; 

for j=k:n
   if k==1
      yvals(:,j) = yvals(:,j) + h*(A*yvals(:,j));
   elseif k==2
      yvals(:,j+1) = yvals(:,j) + (h/2)*(A*(3*yvals(:,j) - yvals(:,j-1)));
   elseif k==3
      yvals(:,j+1) = yvals(:,j) + (h/12)*(A*(23*yvals(:,j) - 16*yvals(:,j-1) + 5*yvals(:,j-2)));
   elseif k==4
      yvals(:,j+1) = yvals(:,j) + (h/24)*(A*(55*yvals(:,j) - 59*yvals(:,j-1) + 37*yvals(:,j-2) - 9*yvals(:,j-3)));
   else
      yvals(:,j+1) = yvals(:,j) + (h/720)*(A*(1901*yvals(:,j) - 2774*yvals(:,j-1) + 2616*yvals(:,j-2) - 1274*yvals(:,j-3) + 251*yvals(:,j-4)));
   end
end
yvals = yvals';


  function [tvals,yvals] = AFixedAM(A,t0,y0,h,k,n)
% [tvals,yvals] = AFixedAM(A,t0,y0,h,k,n)
%
% Let A be a given square matrix.
% Produces an approximate solution to the initial value problem
%
%      y'(t) = Ay(t)     y(t0) = y0
%
% using a strategy that is based upon a k-th order
% Adams-Moulton method. Stepsize is fixed.
%
% A  = square matrix
% t0 = initial time.
% y0 = initial condition vector.
% h  = stepsize.
% k  = order of method. (1<=k<=5).
% n  = number of steps to be taken,
%
% tvals(j) = t0 + (j-1)h, j=1:n+1
% yvals(:j) = approximate solution at t = tvals(j), j=1:n+1

tvals = linspace(t0,t0+h*n,n+1);
yvals = zeros(n+1,length(y0));
[tvals(1:k),yvals(1:k,:)] = ABStart('MyF931',t0,y0,h,k);

yvals = yvals';

I = eye(size(A));
if k==1
   [L,U,P] = lu(I - h*A);
elseif k==2
   [L,U,P] = lu(I - (h/2)*A);
elseif k==3
   [L,U,P] = lu(I - (5*h/12)*A);
elseif k==4
   [L,U,P] = lu(I - (9*h/24)*A);
else
   [L,U,P] = lu(I - (251*h/720)*A);
end
for j=k:n
   if k==1
      rhs = yvals(:,j);        
   elseif k==2
      rhs = yvals(:,j) + (h/2)  *(A*yvals(:,j));
   elseif k==3
      rhs = yvals(:,j) + (h/12) *(A*(8*yvals(:,j) - yvals(:,j-1)));
   elseif k==4
      rhs = yvals(:,j) + (h/24) *(A*(19*yvals(:,j) - 5*yvals(:,j-1) + yvals(:,j-2)));
   else
      rhs = yvals(:,j) + (h/720)*(A*(646*yvals(:,j) - 264*yvals(:,j-1) + 106*yvals(:,j-2) -19*yvals(:,j-3)));
   end
   z = LtriSol(L,P*rhs);
   yvals(:,j+1) = UTriSol(U,z);
end

yvals = yvals';



function yp = MyF931(t,y)
global A
yp = A*y;

P9.3.2

Not Available.

P9.3.3

Not Available.