% scripts and functions for MAT 331, Tuesday, Sept 19, 2017 %------------------------------------------------- % script compare a function to it Taylor series N=20; f=@(x) exp(x); a=-2; b=2; series=1./factorial(0:N); %f =@(x) sin(x); a=-2*pi; b=2*pi; %series= imag(i.^[0:N])./factorial(0:N); %f = @(x) 1./(1-x); a=-1;b=1; %series = ones(1,20); figure; hold on; grid on; title('function and taylor approximations'); x=[a:.01:b]; for k=1:10 plot(x,polyval(series(k:-1:1),x),'LineWidth',2); end plot(x,f(x),'k','LineWidth',3); %------------------------------------------------- % function to compute polynomial interpolation % input is a vector x of distinct real numbers % and a equal length vector y. Finds polynomial % so that p(x_k) = y_k function p = poly_inter(x,y) L =length(x); t=[1:L] p=zeros(1,L); % initialize output to zero for k=1:L s=t(t~=k); % remove kth index L_k=poly(x(s)); % make polynomial with roots at x_j, j~=k L_k = L_k / polyval(L_k,x(k)); % normalize so L_k(x(k))=1 p = p + y(k) *L_k; end % for k return; figure; hold on; grid on; title('Polynomial interpolation'); scatter(x,y,'r'); a=min(x); b=max(x); delta=(b-a)/100; t=[a:delta:b]; plot(t,polyval(p,t),'b'); return; %------------------------------------------------- % function to compute polynomial interpolation % of a given function at evenly points % so that p(x_k) = y_k % try things like exp, sin, abs(x), 1/(1+x^2) function p = even_inter(f,a,b,n) delta=(b-a)/n; x= [a:delta:b]; poly_inter(x,f(x)); return; %------------------------------------------------- % function to compute polynomial interpolation % of a given function at the Chebyshev points function p = cheb_inter(f,a,b,n) x=cos(pi * [0:n]/n) x= a+ (b-a)*(x+1)/2; poly_inter(x,f(x)); return; %------------------------------------------------- % function to use polynomial interpolation to integrate % interpolation points are evenly spaced function out = even_integrate(f,a,b,n) % x=cos(pi * [0:n]/n); % define Chebushev points on [-1,1] % x= a+ (b-a)*(x+1)/2; % rescale x to fit [a,b] delta =(b-a)/n; x=[a:delta:b]; L=length(x); t=[1:L]; y=f(x); for k=1:L s=find(t~=k); % remove kth index from t L_k=poly(x(s)); % make polynomial with roots at x_j, j~=k L_k = L_k / polyval(L_k,x(k)); % normalize so L_k(x(k))=1 Q_k=polyint(L_k); w(k) = diff(polyval(Q_k,[a,b])); end % for k out = sum(w.*f(x)); % now plot the functions (not needed for integral) p=poly_inter(x,y) figure; hold on; grid; title('Plots of function and polynomial with interpolation points'); t=[a:(b-a)/100:b]; plot(t,polyval(p,t),'r','LineWidth',3); % plot poly plot(t,f(t),'--k','LineWidth',2); % plot function scatter(x,polyval(p,x),'k','filled') % plot interpolation poi9nts return %------------------------------------------------- % plot function and its polynomial interpolation function out = plot_inter(f,x) a=min(x); b=max(x); delta=(b-a)/100; p=poly_inter(x,f(x)) figure; hold on; grid; title('Plots of function and polynomial at Chebyshev points'); t=[a:(b-a)/100:b]; plot(t,polyval(p,t),'r','LineWidth',3); % plot poly plot(t,f(t),'--k','LineWidth',2); % plot function scatter(x,y,'k','filled') % plot interpolation poi9nts return %------------------------------------------------- % function to use polynomial interpolation to integrate % interpolation points are Chebyshev points function out = cheb_integrate(f,a,b,n) x=cos(pi * [0:n]/n); % define Chebushev points on [-1,1] x= a+ (b-a)*(x+1)/2; % rescale x to fit [a,b] % delta =(b-a)/n; % used for evenly spaced points % x=[a:delta:b]; L=length(x); t=[1:L]; y=f(x); for k=1:L s=find(t~=k); % remove kth index from t L_k=poly(x(s)); % make polynomial with roots at x_j, j~=k L_k = L_k / polyval(L_k,x(k)); % normalize so L_k(x(k))=1 Q_k=polyint(L_k); w(k) = diff(polyval(Q_k,[a,b])); end % for k out = sum(w.*f(x)); % plot function and polynomials (not needed) plot_inter(f,x); return %------------------------------------------------- %% % Chris Bishop % Homework 1, MAT 331 %% % Set initial variables x=1; N=20; %% % Compute Formula (1) in homework for k=1:N formula2(k)=(1+x/k)^k; end %% % Compute Formula 2: for k=1:N y(k)=x^k/factorial(k); end formula1=1+cumsum(y); %% % Nest we plot both fomulas to show how they both approach e. Note that % Formula 1 appoaches the limit faster figure; hold on; plot(formula1) plot(formula2) figure; hold on; plot(abs(formula1 - exp(x))) plot(abs(formula2 - exp(x))) legend('Formula 1', 'Formula 2');