Lecture Handouts and Examples
| CS 322 Home |
Lec 3 | Feb 2 |
Lec 4 | Feb 4 |
Lec 10 | Feb 25 |
Lecture 10 Wednesday, February 25
% Script Lec10 % Looks at second derivative continuity of pwcubic hermite interpolant bigscreen close all for s2 = [-10 -5 -1 -5/6 1 5 10 100] figure ShowSmallSpline(s2) pause(2) end function ShowSmallSpline(s2) % Displays the piecewise cubic hermite interpolant of the data % % (0,1,2), (1,3,s2), (2,-1,1) % % On [0,1] we have the cubic q1(z) = 1 + 2z +(s2-2)z^2(z-1) which satisfies % % q1(0) = 1, q1'(0) = 2, q1(1) = 3, q1'(3) = s2 % % It has second derivative % % q1"(z) = 2(s2-2) + 6(s2-2)z % The cubic q2(z) = 3 + s2(z-1) + (-4 - s2)(z-1)^2 + (9+s2)(z-1)^2(z-2) % % q2(1) = 3, q1'(1) = s2, q2(2)=-1, q2'(2) = 1 % % It has second derivative % % q2"(z) = 2(-4-s2) + (9+s2)(4(z-1) + 2(z-2)) % % The equation q1"(1) = q2"(1) says % % 2(s2-2) +6(s2-2) = 2(-4-s2) + (9+s2)(-2) % % i.e., s2 = 5/6. When this is the case we have a cubic spline interpolant. close all % Evaluate the two local cubics... zL = linspace(0,1); q1Vals = 1 + 2*zL + (s2-2)*(zL.^2).*(zL-1); zR = linspace(1,2); q2Vals = 3 + s2*(zR-1) + (-4-s2)*(zR-1).^2 + (9+s2)*((zR-1).^2 ).*(zR-2); % Display subplot(2,1,1) plot(zL,q1Vals,zR,q2Vals,[0 1 2],[1 3 -1],'o') legend('q1(z)','q2(z)') title(sprintf('An n=2 piecewise cubic hermite interpolant with s2 = %5.2f',s2),'fontsize',18) % Display their second derivatives subplot(2,1,2) plot(zL,2*(s2-2) + 6*(s2-2)*zL, zR,2*(-4-s2) + (9+s2)*(4*(zR-1)+2*(zR-2))) legend('q1"(z)','q2"(z)') title('It''s second derivative...','fontsize',18)
Lecture 4. Wednesday, February 4
Read about InterpV and HornerV in Chapter 2 of SCMV.
% Lec4A % Plots random cubic interpolants of sin(x) on [0,2pi]. % Uses the Vandermonde method. close all x0 = linspace(0,2*pi,100)'; y0 = sin(x0); for eg=1:30 x = 2*pi*sort(rand(4,1)); y = sin(x); a = InterpV(x,y); pVal = HornerV(a,x0); plot(x0,y0,x0,pVal,'--',x,y,'r*') axis([0 2*pi -2 2]) title('A Random Cubic Interpolant of sin(x) on [0,2\pi]' ) legend('sin(x)','Cubic Interpolant') grid on shg pause end clc grid off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lec4B % Equally spaced interpolants of sin(x) on [0,2pi]. % Uses the Vandermonde method. close all x0 = linspace(0,2*pi,100)'; y0 = sin(x0); for deg=1:20 x = linspace(0,2*pi,deg+1)'; y = sin(x); a = InterpV(x,y); pVal = HornerV(a,x0); subplot(2,1,1) plot(x0,y0,x0,pVal,'--',x,y,'r*') axis([0 2*pi -2 2]) s = sprintf('deg %1d Interpolant of sin(x)',deg); title(s) legend('sin(x)','Polynomial Interpolant') subplot(2,1,2) absErr = abs(pVal-y0); plot(x0,absErr) axis([0 2*pi -inf inf]) s = sprintf('Approx maximum absolute error = %5.1e',max(absErr)); title(s,'fontsize',24) shg pause end clc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lec4c % Equally spaced interpolants of 1/(1+25x^2) on [-1 1]. % Uses the Vandermonde method. close all f = inline('1./(1+25*x.^2)'); x0 = linspace(-1,1,100)'; y0 = feval(f,x0); for deg=1:20 x = linspace(-1,1,deg+1)'; y = feval(f,x); a = InterpV(x,y); pVal = HornerV(a,x0); subplot(2,1,1) plot(x0,y0,x0,pVal,'--',x,y,'r*') s = sprintf('deg %1d Interpolant of 1/(1+25x^{2})',deg); title(s,'fontsize',24) legend('1/(1+25x^{2})','Polynomial Interpolant',0) subplot(2,1,2) absErr = abs(pVal-y0); plot(x0,absErr) s = sprintf('Approx maximum absolute error = %5.1e',max(absErr)); title(s,'fontsize',24) shg pause end clc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Script Lec4d close all ShowInterp('sin',0,2*pi,5) figure g = inline('sin(2*x)'); ShowInterp(g,0,2*pi,5) figure ShowInterp(@sin,0,2*pi,5) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function ShowInterp(fName,a,b,deg) % fName is a string that names a function f defined on % [a,b]. Assume a < b % deg is a positive integer % Displays the equal-spacing, degree d polynomial interpolant % of f on [a,b]. % Display the function over the interval... x = linspace(a,b)'; y = feval(fName,x); plot(x,y); axis([a b -inf inf]) % Display the points of interpolation hold on x0 = linspace(a,b,deg+1)'; y0 = feval(fName,x0); plot(x0,y0,'or') % Generate the interpolant a = InterpV(x0,y0); % Evaluate the interpolant across the interval and display. pVals = HornerV(a,x); plot(x,pVals,'r') hold off shg
Lecture 3. Monday, February 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Script Lec3A x = linspace(-7,7); y = MyExp(x); shg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = MyExp(x) % x is a vector % y has the same length and orientation with % y(i) = exp(x(i)), i = 1:length(x) y = ones(size(x)); % y = vector of all 1's, same size & shape as x termIndex = 1; termNext = x; yNext = y + termNext; while(any(y~=yNext)) % Iterate as long as the addition of the next term makes a difference y = yNext; termIndex = termIndex + 1; termNext = (termNext.*x)/termIndex; yNext = yNext + termNext; % Display the relative error relErr = abs(yNext - exp(x))./exp(x); % Add eps to relErr so don't evaluate log at zero. semilogy(x,relErr+eps,[min(x) max(x)] ,[10^-15 10^-15],'r') title(sprintf('Number of terms = %1d',termIndex)) ylabel('Relative Error') xlabel('x') shg pause end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Script Lec3B % Explores running sums in Taylor expansion for exp(x) close all x = -7; nTerms = 30; sum = zeros(1,nTerms); term = zeros(1,nTerms); sum(1) = 1 + x; term(1) = x; for k=2:nTerms term(k) = x*term(k-1)/k; sum(k) = sum(k-1) + term(k); end semilogy(1:nTerms,abs(sum),1:nTerms,abs(term)) title(sprintf('exp(%5.1f) via taylor series',x)) xlabel('k') legend('1 + term(1) + ... + term(k)',' | term(k) |') shg