Chapter 9 M-Files
Script Files
ShowTraj | Shows family of solutions. |
ShowEuler | Illustrates Euler method. |
ShowFixedEuler | Plots error in fixed step Euler for y'=y, y(0)=1. |
ShowTrunc | Shows effect of truncation error. |
EulerRoundoff | Illustrates Euler in three-digit floating point. |
ShowAB | Illustrates FixedAB. |
ShowPC | Illustrates FixedPC. |
ShowRK | Illustrates FixedRK. |
ShowKepler | Illustrates ode23 and ode45 on a system. |
Function Files
FixedEuler | Fixed step Euler method. |
ABStart | Gets starting values for Adams methods. |
ABStep | Adams-Bashforth step (order <=5). |
FixedAB | Fixed step size Adams-Bashforth. |
AMStep | Adams-Moulton step (order <=5). |
PCStep | AB-AM predictor-corrector Step (order <=5$). |
FixedPC | Fixed stepsize AB-AM predictor-corrector. |
RKStep | Runge-Kutta step (order <=5). |
FixedRK | Fixed step size Runge-Kutta. |
Kepler | For solving two-body IVP. |
f1 | The f function for the model problem y'=y. |
% 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)')
% 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
% 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)))
% 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')
% 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')
% 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])
% 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])
% 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')
% 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')
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
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
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);
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
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);
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);
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
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);
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
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] ;
function yp = f1(t,y) yp = -y;