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.