Chapter 9: Problem Solutions
% 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;
% 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);
% 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;
Not Available
% 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;
% 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
% 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;
% 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));
Not Available
% 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);
Noty Available
% 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;
Not Available.
Not Available.