Chapter 1 M-Files
Script File Index
SineTable | Prints a short table of sine evaluations. |
SinePlot | Displays a sequence of sin(x) plots. |
ExpPlot | Plots exp(x) and an approximation to exp(x). |
TangentPlot | Plots tan(x). |
SineAndCosPlot | Superimposes plots of sin(x) and cos(x). |
Polygons | Displays nine regular polygons, one per window. |
SumOfSines | Displays the sum of four sine functions. |
SumOfSines2 | Displays a pair of sum-of-sine functions. |
UpDown | Sample core exploratory environment. |
RunUpDown | Framework for running UpDown. |
Histograms | Displays the distribution of rand and randn. |
Clouds | Displays 2-dimensional rand and randn. |
Dice | Histogram of 1000 dice rolls. |
Darts | Monte Carlo computation of pi. |
Smooth | Polygon smoothing. |
Stirling | Relative and absolute error in Stirling formula. |
ExpTaylor | Plots relative error in Taylor approximation to exp(x). |
Zoom | Roundoff in the expansion of (x-1)^ 6. |
FpFacts | Examines precision, overflow, and underflow. |
TestMyExp | Examines MyExp1 , MyExp2 , MyExp3 , and MyExp4. |
TestDerivative | Applies Derivative to Sin10. |
Euler | Three-digit arithmetic sum of 1 + 1/2 + ... + 1/n. |
ShowPadeArray | Tests the function PadeArray. |
ShowFonts | Illustrates how to use fonts. |
ShowSymbols | Shows how to generate math symbols. |
ShowGreek | Shows how to generate Greek letters. |
ShowText | Shows how to align with text. |
ShowLineWidth | Shows how vary line width in a plot. |
ShowAxes | Shows how to set tick marks on axes. |
ShowLegend | Shows how to add a legend to a plot. |
ShowColor | Shows how to use built-in colors and user-defined colors. |
Function File Index
MyExpF | For-loop Taylor approximation to exp(x). |
MyExp1 | Vectorized version of MyExpF. |
MyExp2 | Better vectorized version of MyExpF. |
MyExpW | While-loop Taylor approximation to exp(x). |
MyExp3 | Vectorized version of MyExpW. |
MyExp4 | Better vectorized version of MyExpW. |
Derivative | Numerical differentiation. |
Sin10 | sin(10x). |
Represent | Sets up 3-digit arithmetic representation. |
Convert | Converts 3-digit representation to float. |
Float | Simulates 3-digit arithmetic. |
Pretty | Pretty prints a 3-digit representation. |
PadeArray | Builds a cell array of Pade coefficients. |
PadeCoeff | Generates Pade approximant coefficients |
% Script File: SineTable % Prints a short table of sine evaluations. clc n = 21; x = linspace(0,1,n); y = sin(2*pi*x); disp(' ') disp(' k x(k) sin(x(k))') disp('------------------------') for k=1:21 degrees = (k-1)*360/(n-1); disp(sprintf(' %2.0f %3.0f %6.3f ',k,degrees,y(k))); end disp( ' '); disp('x(k) is given in degrees.') disp(sprintf('One Degree = %5.3e Radians',pi/180))
% Script File: SinePlot % Displays increasingly smooth plots of sin(2*pi*x). close all for n = [4 8 12 16 20 50 100 200 400] x = linspace(0,1,n); y = sin(2*pi*x); plot(x,y) title(sprintf('Plot of sin(2*pi*x) based upon n = %3.0f points.',n)) pause(1) end
% Script File: ExpPlot % Examines the function % % f(x) = ((1 + x/24)/(1 - x/24 + x^2/384))^8 % % as an approximation to exp(z) across [0,1]. close all x = linspace(0,1,200); num = 1 + x/24; denom = 1 - x/12 + (x/384).*x; quot = num./denom; y = quot.^8; plot(x,y,x,exp(x))
% Script File: TangentPlot % Plots the function tan(x), -pi/2 <= x <= 9pi/2 close all ymax = 10; x = linspace(-pi/2,pi/2,40); y = tan(x); plot(x,y) axis([-pi/2 9*pi/2 -ymax ymax]) title('The Tangent Function') xlabel('x') ylabel('tan(x)') hold on for k=1:4 xnew = x+ k*pi; plot(xnew,y); end hold off
% Script File: SineAndCosPlot % Plots the functions sin(2*pi*x) and cos(2*pi*x) across [0,1] % and marks their intersection. close all x = linspace(0,1,200); y1 = sin(2*pi*x); y2 = cos(2*pi*x); plot(x,y1,x,y2,'--',[1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')
% Script File: Polygons % Plots selected regular polygons. close all theta = linspace(0,2*pi,361); c = cos(theta); s = sin(theta); k=0; for sides = [3 4 5 6 8 10 12 18 24] stride = 360/sides; k=k+1; subplot(3,3,k) plot(c(1:stride:361),s(1:stride:361)) axis([-1.2 1.2 -1.2 1.2]) axis equal end
% Script File: SumOfSines % Plots f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) + 5sin(4x) % across the interval [-10,10]. close all x = linspace(-10,10,200)'; A = [sin(x) sin(2*x) sin(3*x) sin(4*x)]; y = A*[2;3;7;5]; plot(x,y) title('f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) + 5sin(4x)')
% Script File: SumOfSines2 % Plots the functions % f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) + 5sin(4x) % g(x) = 8sin(x) + 2sin(2x) + 6sin(3x) + 9sin(4x) % across the interval [-10,10]. close all n = 200; x = linspace(-10,10,n)'; A = [sin(x) sin(2*x) sin(3*x) sin(4*x)]; y = A*[2 8;3 2;7 6;5 9]; plot(x,y)
% Script File: UpDown % Generates a column vector x(1:n) of positive integers % where x(1) is solicited and % % x(k+1) = x(k)/2 if x(k) is even. % x(k+1) = 3x(k)+1 if x(k) is odd. % % The value of n is either 500 or the first index with the % property that x(n) = 1, whichever comes first. x = zeros(500,1); x(1) = input('Enter initial positive integer:'); k = 1; while ((x(k) ~= 1) & (k < 500)) if rem(x(k),2) == 0 x(k+1) = x(k)/2; else x(k+1) = 3*x(k)+1; end k = k+1; end n = k; x = x(1:n); clc disp(sprintf('x(1:%1.0f) = \n',n)) disp(sprintf('%-8.0f',x)) [xmax,imax] = max(x); disp(sprintf('\n x(%1.0f) = %1.0f is the max.',imax,xmax)) density = sum(x<=x(1))/x(1); disp(sprintf(' The density is %5.3f.',density)) close all figure plot(x) title(sprintf('x(1) = %1.0f, n = %1.0f',x(1),n)); figure plot(-sort(-x)) title('Sequence values sorted.') I = find(rem(x(1:n-1),2)==1); if length(I)>1 figure plot((1:n),zeros(1,n),I+1,x(I+1),I+1,x(I+1),'*') title('Local Maxima') end
% Script File: RunUpDown % Environment for studying the up/down sequence. % Stores selected results in file UpDownOutput. while(input('Another Example? (1=yes, 0=no)')) diary UpDownOutput UpDown diary off disp(' ') if (input('Keep Output? (1=yes, 0=no)')~=1) delete UpDownOutput end end
% Script File: Histograms % Histograms of rand(1000,1) and randn(1000,1). close all subplot(2,1,1) x = rand(1000,1); hist(x,30) axis([-1 2 0 60]) title('Distribution of Values in rand(1000,1)') xlabel(sprintf('Mean = %5.3f. Median = %5.3f.',mean(x),median(x))) subplot(2,1,2) x = randn(1000,1); hist(x,linspace(-2.9,2.9,100)); title('Distribution of Values in randn(1000,1)') xlabel(sprintf('Mean = %5.3f. Standard Deviation = %5.3f',mean(x),std(x)))
% Script File: Clouds % 2-dimensional pictures of the uniform and normal distributions. close all Points = rand(1000,2); subplot(1,2,1) plot(Points(:,1),Points(:,2),'.') title('Uniform Distribution.') axis([0 1 0 1]) axis square Points = randn(1000,2); subplot(1,2,2) plot(Points(:,1),Points(:,2),'.') title('Normal Distribution.') axis([-3 3 -3 3]) axis square
% Script File: Dice % Simulates 1000 rollings of a pair of dice. close all First = 1 + floor(6*rand(1000,1)); Second = 1 + floor(6*rand(1000,1)); Throws = First + Second; hist(Throws, linspace(2,12,11)); title('Outcome of 1000 Dice Rolls.')
% Script File: Darts % Estimates pi using random dart throws. close all rand('seed',.123456) NumberInside = 0; PiEstimate = zeros(500,1); for k=1:500 x = -1+2*rand(100,1); y = -1+2*rand(100,1); NumberInside = NumberInside + sum(x.^2 + y.^2 <= 1); PiEstimate(k) = (NumberInside/(k*100))*4; end plot(PiEstimate) title(sprintf('Monte Carlo Estimate of Pi = %5.3f',PiEstimate(500))); xlabel('Hundreds of Trials')
% Script File: Smooth % Solicits n, draws an n-gon, and then smooths it. close all n = input('Enter the number of edges:'); figure axis([0 1 0 1]) axis square hold on x = zeros(n,1); y = zeros(n,1); for k=1:n title(sprintf('Click in %2.0f more points.',n-k+1)) [x(k) y(k)] = ginput(1); plot(x(1:k),y(1:k), x(1:k),y(1:k),'*') end x0 = x; y0 = y; x = [x;x(1)]; y = [y;y(1)]; plot(x,y,x,y,'*') title('The Original Polygon') hold off k=0; xlabel('Click inside window to smooth, outside window to quit.') [a,b] = ginput(1); v = axis; while (v(1)<=a) & (a<=v(2)) & (v(3)<=b) & (b<=v(4)); k = k+1; x = [(x(1:n)+x(2:n+1))/2;(x(1)+x(2))/2]; y = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2]; m = max(abs([x;y])); x = x/m; y = y/m; figure plot(x,y,x,y,'*') axis square title(sprintf('Number of Smoothings = %1.0f',k)) xlabel('Click inside window to smooth, outside window to quit.') v = axis; [a,b] = ginput(1); end
% Script File: Stirling % Prints a table showing error in Stirling's formula for n! close all clc disp( ' ') disp(' Stirling Absolute Relative') disp(' n n! Approximation Error Error') disp('----------------------------------------------------------------') e=exp(1); nfact=1; for n=1:13 nfact = n*nfact; s = sqrt(2*pi*n)*((n/e)^n); abserror = abs(nfact - s); relerror = abserror/nfact; s1 = sprintf(' %2.0f %10.0f %13.2f',n,nfact,s); s2 = sprintf(' %13.2f %5.2e',abserror,relerror); disp([s1 s2]) end
% Script File: ExpTaylor % Plots, as a function of n, the relative error in the % Taylor approximation % % 1 + x + x^2/2! +...+ x^n/n! % % to exp(x). close all nTerms = 50; for x=[10 5 1 -1 -5 -10] figure f = exp(x)*ones(nTerms,1); s = 1; term = 1; for k=1:50 term = x.*term/k; s = s+ term; err(k) = abs(f(k) - s); end relerr = err/exp(x); semilogy(1:nTerms,relerr) ylabel('Relative Error in Partial Sum.') xlabel('Order of Partial Sum.') title(sprintf('x = %5.2f',x)) end
% Script File: Zoom % Plots (x-1)^6 near x=1 with increasingly refined scale. % Evaluation via x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x +1 % leads to severe cancellation. close all k=0; n=100; for delta = [.1 .01 .008 .007 .005 .003 ] x = linspace(1-delta,1+delta,n)'; y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 - 6*x + ones(n,1); k = k+1; subplot(2,3,k) plot(x,y,x,zeros(1,n)) axis([1-delta 1+delta -max(abs(y)) max(abs(y))]) end
% Script File: FpFacts % Prints some facts about the underlying floating point system. close all clc % p = smallest positive integer so 1+1/2^p = 1. x=1; p=0; y=1; z=x+y; while x~=z y=y/2; p=p+1; z=x+y; end disp(sprintf('p = %2.0f is the smallest positive integer so 1+1/2^p = 1.',p)) % q = smallest positive integer so 1/2^q = 0. x=1; q=0; while x>0 x=x/2; q=q+1; end; disp(sprintf('q = %2.0f is the smallest positive integer so 1/2^q = 0.',q)) % r = smallest positive integer so 2^r = inf. x=1; r=0; while x~=inf x=2*x; r=r+1; end disp(sprintf('r = %2.0f is the smallest positive integer so 2^r = inf.',r)) % The Kahan example. In exact arithmetic f/e = 0/0: disp(' ') echo on h = 1./2.; x = 2./3. - h; y = 3./5. - h; e = (x+x+x) - h; f = (y+y+y+y+y)-h; z = f/e; echo off disp(sprintf('\nz = %10.5f',z))
% Script File: TestMyExp % Checks MyExpF, MyExp1, MyExp2, MyExp3, and MyExp4. close all m = 50; x = linspace(-1,1,m); exact = exp(x); figure k=0; y = zeros(1,m); for n = [4 8 16 20] for i=1:m y(i) = MyExpF(x(i),n); end RelErr = abs(exact - y)./exact; k=k+1; subplot(2,2,k) plot(x,RelErr) title(sprintf('n = %2.0f',n)) end clc nRepeat = 1; % Note: The following fragment may take a long time to execute % with this value of nRepeat. disp('T1 = MyExp1 benchmark.') disp('T2 = MyExp2 benchmark.') disp(' ') disp(' Length(x) T2/T1') disp('-----------------------') for L = 1000:100:2000 xL = linspace(-1,1,L); tic for k=1:nRepeat y = MyExp1(xL); end T1 = toc; tic for k=1:nRepeat y = MyExp2(xL); end T2 = toc; disp(sprintf('%6.0f %13.6f ',L,T2/T1)) end disp(' ') disp(' ') disp(' Length(x) MyExp3(x) MyExp4(x) ') disp(' Flops Flops') disp('---------------------------------------') for L = 50:50:300 xL = linspace(-1,1,L); flops(0); y = MyExp3(xL); f3 = flops; flops(0); y = MyExp4(xL); f4 = flops; disp(sprintf('%6.0f %13.0f %13.0f ',L,f3,f4)) end
% Script File: TestDerivative % Numerical differentiation of f(x) = sin(10x) close all a = 10; M2=100; % Look at errors across [0,1]. % err1 = error with no derivative information. % err2 = error with 2nd derivative information. m = 100; x=linspace(0,1,m); for i=1:m exactDer = a*cos(a*x(i)); err1(i) = abs(exactDer - Derivative('sin10',x(i))); err2(i) = abs(exactDer - Derivative('sin10',x(i),eps,M2)); end % Plot err1(1:m) figure plot(x,err1); title('Derivative(''sin10'',a)') xlabel('a') ylabel('Der Error') % Plot err2(1:m) figure plot(x,err2); title('Derivative(''sin10'',a,eps,100)') xlabel('a') ylabel('Der Error')
% Script File: Euler % Sums the series 1 + 1/2 + 1/3 + .. in 3-digit floating point arithmetic. % Terminates when the addition of the next term does not change % the value of the running sum. oldsum = represent(0); one = represent(1); sum = one; k = 1; while convert(sum) ~= convert(oldsum) k = k+1; kay = represent(k); term = float(one,kay,'/'); oldsum = sum; sum = float(sum,term,'+'); end clc disp(['The sum for ' num2str(k) ' or more terms is ' pretty(sum)])
% Script File: ShowPadeArray % Illustrates the function PadeArray that creates a cell array. clc P = PadeArray(5,5); for i=1:5 for j=1:5 disp(sprintf('\n (%1d,%1d) Pade Coefficients:\n',i-1,j-1)) disp([' ' sprintf('%7.4f ',P{i,j}.num)]) disp([' ' sprintf('%7.4f ',P{i,j}.den)]) end end
% Script File: ShowFonts close all HA = 'HorizontalAlign'; fonts = {'Times-Roman' 'Helvetica' 'AvantGarde' 'Bookman' 'Palatino'... 'ZapfChancery' 'Courier' 'NewCenturySchlbk' 'Helvetica-Narrow'}; for k=1:9 figure plot([-10 100 100 -10 -10],[0 0 60 60 0]) axis([-10 100 0 60]) axis off v=38; F = fonts{k}; text(45,55,F,'FontName',F,'FontSize',24,HA,'center') text(10,47,'Plain','FontName',F,'FontSize',22,HA,'center') text(45,47,'Bold','FontName',F,'Fontweight','bold','FontSize',22,HA,'center') text(82,47,'Oblique','FontName',F,'FontAngle','oblique','FontSize',22,HA,'center') for size=[22 18 14 12 11 10 9] text(10,v,'Matlab','FontName',F,'FontSize',size,HA,'center') text(45,v,'Matlab','FontName',F,'FontSize',size,HA,'center','FontWeight','bold') text(82,v,'Matlab','FontName',F,'FontSize',size,HA,'center','FontAngle','oblique') v = v-6; end end
% Script File: ShowSymbols close all plot([0 12 12 0 0],[0 0 12 12 0]) axis off text(6,10.5,'Math Symbols','FontSize',18,'HorizontalAlign','center') x = 1; x1 = x+.7; y = 4.6; y1 = y+.7; z = 9; z1 = z+.7; n = 12; text(y,9,'\leftarrow','Fontsize',n); text(y,8,'\rightarrow','Fontsize',n); text(y,7,'\uparrow','Fontsize',n) text(y,6,'\downarrow','Fontsize',n) text(y,5,'\Leftarrow','FontSize',n) text(y,4,'\Rightarrow','FontSize',n) text(y,3,'\Leftrightarrow','FontSize',n) text(y,2,'\partial','FontSize',n) text(x,9,'\neq','FontSize',n) text(x,8,'\geq','FontSize',n) text(x,7,'\approx','FontSize',n) text(x,6,'\equiv','Fontsize',n) text(x,5,'\cong','Fontsize',n) text(x,4,'\pm','Fontsize',n) text(x,3,'\nabla','FontSize',n) text(x,2,'\angle','FontSize',n) text(z,9,'\in','FontSize',n) text(z,8,'\subset','Fontsize',n) text(z,7,'\cup','Fontsize',n) text(z,6,'\cap','Fontsize',n) text(z,5,'\perp','FontSize',n) text(z,4,'\infty','FontSize',n) text(z,3,'\int','Fontsize',n) text(z,2,'\times','FontSize',n) text(y1,9,'\\leftarrow','Fontsize',n); text(y1,8,'\\rightarrow','Fontsize',n); text(y1,7,'\\uparrow','Fontsize',n) text(y1,6,'\\downarrow','Fontsize',n) text(y1,5,'\\Leftarrow','FontSize',n) text(y1,4,'\\Rightarrow','FontSize',n) text(y1,3,'\\Leftrightarrow','FontSize',n) text(y1,2,'\\partial','FontSize',n) text(x1,9,'\\neq','FontSize',n) text(x1,8,'\\geq','FontSize',n) text(x1,7,'\\approx','FontSize',n) text(x1,6,'\\equiv','Fontsize',n) text(x1,5,'\\cong','Fontsize',n) text(x1,4,'\\pm','Fontsize',n) text(x1,3,'\\nabla','FontSize',n) text(x1,2,'\\angle','FontSize',n) text(z1,9,'\\in','FontSize',n) text(z1,8,'\\subset','Fontsize',n) text(z1,7,'\\cup','Fontsize',n) text(z1,6,'\\cap','Fontsize',n) text(z1,5,'\\perp','FontSize',n) text(z1,4,'\\infty','FontSize',n) text(z1,3,'\\int','Fontsize',n) text(z1,2,'\\times','FontSize',n)
% Script File: ShowGreek close all plot([-1 12 12 -1 -1],[-1 -1 12 12 -1]) axis off text(4,10,'Greek Symbols','FontSize',18) x = 0; x1 = x+.7; y = 4; y1 = y+.7; z = 8; z1 = z+.7; n = 12; text(x,8,'\alpha','Fontsize',n); text(x,7,'\beta','Fontsize',n); text(x,6,'\gamma','Fontsize',n) text(x,5,'\delta','Fontsize',n) text(x,4,'\epsilon','FontSize',n) text(x,3,'\kappa','FontSize',n) text(x,2,'\lambda','Fontsize',n) text(x,1,'\mu','Fontsize',n) text(x,0,'\nu','Fontsize',n) text(x1,8,'\\alpha','Fontsize',n); text(x1,7,'\\beta','Fontsize',n); text(x1,6,'\\gamma','Fontsize',n) text(x1,5,'\\delta','Fontsize',n) text(x1,4,'\\epsilon','FontSize',n) text(x1,3,'\\kappa','FontSize',n) text(x1,2,'\\lambda','Fontsize',n) text(x1,1,'\\mu','Fontsize',n) text(x1,0,'\\nu','Fontsize',n) text(y,8,'\omega','Fontsize',n) text(y,7,'\phi','Fontsize',n) text(y,6,'\pi','FontSize',n) text(y,5,'\chi','Fontsize',n) text(y,4,'\psi','Fontsize',n) text(y,3,'\rho','FontSize',n) text(y,2,'\sigma','Fontsize',n) text(y,1,'\tau','FontSize',n) text(y,0,'\upsilon','FontSize',n) text(y1,8,'\\omega','Fontsize',n) text(y1,7,'\\phi','Fontsize',n) text(y1,6,'\\pi','FontSize',n) text(y1,5,'\\chi','Fontsize',n) text(y1,4,'\\psi','Fontsize',n) text(y1,3,'\\rho','FontSize',n) text(y1,2,'\\sigma','Fontsize',n) text(y1,1,'\\tau','FontSize',n) text(y1,0,'\\upsilon','FontSize',n) text(z,8,'\Sigma','FontSize',n) text(z,7,'\Pi','FontSize',n) text(z,6,'\Lambda','FontSize',n) text(z,5,'\Omega','FontSize',n) text(z,4,'\Gamma','FontSize',n) text(z1,8,'\\Sigma','FontSize',n) text(z1,7,'\\Pi','FontSize',n) text(z1,6,'\\Lambda','FontSize',n) text(z1,5,'\\Omega','FontSize',n) text(z1,4,'\\Gamma','FontSize',n)
% Script File: ShowText close all r = 1; t = pi/6 + linspace(0,2*pi,7); x = r*cos(t); y = r*sin(t); plot(x,y); axis equal off HA = 'HorizontalAlignment'; VA = 'VerticalAlignment'; text(x(1),y(1),'\leftarrow {\itP}_{1}', HA,'left','FontSize',14) text(x(2),y(2),'\downarrow', HA,'center',VA,'baseline','FontSize',14) text(x(2),y(2),'{ \itP}_{2}', HA,'left',VA,'bottom','FontSize',14) text(x(3),y(3),'{\itP}_{3} \rightarrow', HA,'right','FontSize',14) text(x(4),y(4),'{\itP}_{4} \rightarrow', HA,'right','FontSize',14) text(x(5),y(5),'\uparrow', HA,'center',VA,'top','FontSize',14) text(x(5),y(5),'{\itP}_{5} ', HA,'right',VA,'top','FontSize',14) text(x(6),y(6),'\leftarrow {\itP}_{6}', HA,'left','FontSize',14) text(0,1.4*r,'A Labeled Hexagon^{1}',HA,'center','FontSize',14) text(0,-1.4*r,'^{1} A hexagon has six sides.',HA,'center','FontSize',10)
% Script File: ShowLineWidth close all plot([0 14 14 0 0],[0 0 14 14 0]) text(7,13,'LineWidth','FontSize',18,'HorizontalAlign','center') axis off hold on for width=0:10 h = plot([3 12],[11-width 11-width]); if width>0 set(h,'LineWidth',width) text(1,11-width,sprintf('%3d',width),'FontSize',14) else text(1,11-width,'default','FontSize',14)
% Script File: ShowAxes F = 'Times-Roman'; n = 12; close all t = linspace(0,2*pi); c = cos(t); s = sin(t); plot(c,s) axis([-1.3 1.3,-1.3 1.3]) axis equal title('The Circle ({\itx-a})^{2} + ({\ity-b})^{2} = {\itr}^{2}','FontName',F,'FontSize',n) xlabel('x','FontName',F,'FontSize',n) ylabel('y','FontName',F,'FontSize',n) set(gca,'XTick',[-.5 0 .5]) set(gca,'YTick',[-.5 0 .5]) grid on
% Script File: ShowLegend close all t = linspace(0,2); axis([0 2 -1.5 1.5]) y1 = sin(t*pi); y2 = cos(t*pi); plot(t,y1,t,y2,'--',[0 .5 1 1.5 2],[0 0 0 0 0],'o') set(gca,'XTick',[]) set(gca,'YTick',[0]) grid on legend('sin(\pi t)','cos(\pi t)','roots',0)
% Script File: ShowColor close all plot([-1 16 16 -1 -1],[-3 -3 18 18 -3]) axis off text(0.5,16.5,'Built-In Colors','FontSize',14) text(9.5,16.5,'A Gradient','FontSize',14) % The predefined colors. x = [3 6 6 3 3]; y = [-.5 -.5 .5 .5 -.5]; hold on fill(x,y,'c') text(0,0,'cyan','FontSize',12) fill(x,y+2,'m') text(0,2,'magenta','FontSize',12) fill(x,y+4,'y') text(0,4,'yellow','FontSize',12) fill(x,y+6,'r') text(0,6,'red','FontSize',12) fill(x,y+8,'g') text(0,8,'green','FontSize',12) fill(x,y+10,'b') text(0,10,'blue','FontSize',12) fill(x,y+12,'k') text(0,12,'black','FontSize',12) fill(x,y+14,'w') text(0,14,'white','FontSize',12) % Making your own colors. a = -.5; b = 14.5; L = 11; R = 15; n = 100; k=1; for f = [.4 .6 .7 .8 .85 .90 .95 1] color = [f 1 1]; text(x(1)+5,y(1)+(2*(k-1)),sprintf('[ %4.2f , %1d , %1d ]',f,1,1),'VerticalAlign','bottom') fill(x+9,y+(2*(k-1)),color) k = k+1; end hold off
function y = MyExpF(x,n) % y = MyExpF(x,n) % x is a scalar, n is a positive integer % and y = n-th order Taylor approximation to exp(x). y = 1; term = 1; for k = 1:n term = x*term/k; y = y + term; end
function y = MyExp1(x) % y = MyExp1(x) % x is a column n-vector and y is a column n-vector with the property that % y(i) is a Taylor approximation to exp(x(i)) for i=1:n. nTerms = 17; n = length(x); y = ones(n,1); for i=1:n y(i) = MyExpF(x(i),nTerms); end
function y = MyExp2(x) % y = MyExp2(x) % x is an n-vector and y is an n-vector with with the same shape % and the property that y(i) is a Taylor approximation to % exp(x(i)) for i=1:n. y = ones(size(x)); nTerms = 17; term = ones(size(x)); for k=1:nTerms term = x.*term/k; y = y + term; end
function y = MyExpW(x) % y = MyExpW(x) % x is a scalar and y is a Taylor approximation to exp(x). y = 0; term = 1; k=0; while abs(term) > eps*abs(y) k = k + 1; y = y + term; term = x*term/k; end
function y = MyExp3(x) % y = MyExp3(x) % x is a column n-vector and y is a column n-vector with the property that % y(i) is a Taylor approximation to exp(x(i)) for i=1:n. n = length(x); y = ones(n,1); for i=1:n y(i) = MyExpW(x(i)); end
function y = MyExp4(x) % y = MyExp4(x) % x is an n-vector and y is an n-vector with the same shape and the % property that y(i) is a Taylor approximation to exp(x(i)) for i=1:n. y = zeros(size(x)); term = ones(size(x)); k = 0; while any(abs(term) > eps*abs(y)) y = y + term; k = k+1; term = x.*term/k; end
function [d,err] = Derivative(fname,a,delta,M2) % d = Derivative(fname,a,delta,M2); % fname a string that names a function f(x) whose derivative % at x = a is sought. delta is the absolute error associated with % an f-evaluation and M2 is an estimate of the second derivative % magnitude near a. d is an approximation to f'(a) and err is an estimate % of its absolute error. % % Usage: % [d,err] = Derivative(fname,a) % [d,err] = Derivative(fname,a,delta) % [d,err] = Derivative(fname,a,delta,M2) if nargin &= 3 % No derivative bound supplied, so assume the % second derivative bound is 1. M2 = 1; end if nargin == 2 % No function evaluation error supplied, so % set delta to eps. delta = eps; end % Compute optimum h and divided difference hopt = 2*sqrt(delta/M2); d = (feval(fname,a+hopt) - feval(fname,a))/hopt; err = 2*sqrt(delta*M2);
function y = sin10(x) y = sin(10*x);
function f = Represent(x) % f = Represent(x) % Yields a 3-digit mantissa floating point representation of f: % % f.mSignBit mantissa sign bit (0 if x>=0, 1 otherwise) % f.m mantissa (= f.m(1)/10 + f.m(2)/100 + f.m(3)/1000) % f.eSignBit the exponent sign bit (0 if exponent nonnegative, 1 otherwise) % f.e the exponent (-9<=f.e<=9) % % If x is out side of [-.999*10^9,.999*10^9], f.m is set to inf. % If x is in the range [-.100*10^-9,.100*10^-9] f is the representation of zero % in which both sign bits are 0, e is zero, and m = [0 0 0]. f = struct('mSignBit',0,'m',[0 0 0],'eSignBit',0,'e',0); if abs(x)<.100*10^-9 % Underflow. Return 0 return end if x>.999*10^9 % Overflow. Return inf f.m = inf; return end if x<-.999*10^9 % Overflow. Return -inf f.mSignBit = 1; f.m = inf; return end % Set the mantissa sign bit if x>0 f.mSignBit = 0; else f.mSignBit = 1; end % Determine m and e so .1 <= m < 1 and abs(x) = m*10^e m = abs(x); e = 0; while m >= 1, m = m/10; e = e+1; end while m < 1/10,m = 10*m; e = e-1; end % Determine nearest integer to 1000m z = round(1000*m); % Set the mantissa f.m(1) = floor(z/100); f.m(2) = floor((z - f.m(1)*100)/10); f.m(3) = z - f.m(1)*100 - f.m(2)*10; % Set the exponent and its sign bit. if e>=0 f.eSignBit = 0; f.e = e; else f.eSignBit = 1; f.e = -e; end
function x = Convert(f) % x = Convert(f) % f is a is a representation of a 3-digit floating point number. (For details % type help represent. x is the value of v. % Overflow situations if (f.m == inf) & (f.mSignBit==0) x = inf; return end if (f.m == inf) & (f.mSignBit==1) x = -inf; return end % Mantissa value mValue = (100*f.m(1) + 10*f.m(2) + f.m(3))/1000; if f.mSignBit==1 mValue = -mValue; end % Exponent value eValue = f.e; if f.eSignBit==1 eValue = -eValue; end x = mValue * 10^eValue;
function z = Float(x,y,op) % z = Float(x,y,op) % x and y are representations of a 3-digit floating point number. (For details % type help represent. % op is one of the strings '+', '-', '*', or '/'. % z is the 3-digit floating point representation of x op y. sx = num2str(convert(x)); sy = num2str(convert(y)); z = represent(eval(['(' sx ')' op '(' sy ')' ]));
function s = Pretty(f) % s = Pretty(f) % f is a is a representation of a 3-digit floating point number. (For details % type help represent. % s is a string so that disp(s) "pretty prints" the value of f. if (f.m == inf) & (f.mSignBit==0) s = 'inf'; elseif (f.m==inf) & (f.mSignBit==1) s = '-inf'; else % Convert the mantissa. m = ['.' num2str(f.m(1)) num2str(f.m(2)) num2str(f.m(3))]; if f.mSignBit == 1 m = ['-' m]; end % Convert the exponent. e = num2str(f.e); if f.eSignBit == 0 e = ['+' num2str(f.e)]; else e = ['-' num2str(f.e)]; end s = [m ' x 10^' e ]; end
function P = PadeArray(m,n) % P = PadeArray(m,n) % m and n are nonnegative integers. % P is an (m+1)-by-(n+1) cell array. % % P{i,j} represents the (i-1,j-1) Pade approximation N(x)/D(x) to exp(x). P = cell(m+1,n+1); for i=1:m+1 for j=1:n+1 P{i,j} = PadeCoeff(i-1,j-1); end end
function R = PadeCoeff(p,q) % R = PadeCoeff(p,q) % p and q are nonnegative integers and R is a representation of the % (p,q)-Pade approximation N(x)/D(x) to exp(x): % % R.num is a row (p+1)-vector whose entries are the coefficients of the % p-degree numerator polynomial N(x). % % R.den is a row (q+1)-vector whose entries are the coefficients of the % q-degree denominator polynomial D(x). % % Thus, % R.num(1) + R.num(2)x + R.num(3)x^2 % ------------------------------------ % R.den(1) + R.den(2)x % % is the (2,1) Pade approximation. a(1) = 1; for k=1:p, a(k+1) = a(k)*(p-k+1)/(k*(p+q-k+1)); end b(1) = 1; for k=1:q, b(k+1) = -b(k)*(q-k+1)/(k*(p+q-k+1)); end R = struct('num',a,'den',b);