Math 104A Homework #2 Instructor: Lihui Chai General Instructions: Please write your homework papers neatly. You need to turn in both full printouts of your codes and the appropriate runs you made. Write your own code, individually. Do not copy codes! 1. (a) Write the Lagrangian form of the interpolating polynomial P2 (x) corresponding to the data in the table below: xj 0 1 3 f (xj ) 1 1 -5 (b) Use P2 (x) you obtained in (a) to approximate f (2). 2. We proved in class that kf Pn k1 , Pn k1 (1 + n ) kf (1) where Pn is the interpolating polynomial of f at the nodes x0 , ..., xn , Pn is the best approximation of f , in the maximum (infinity) norm, by a polynomial of degree at most n, and n = n X (n) lj j=0 (n) is the Lebesgue constant (here the lj , (2) 1 are the elementary Lagrange polynomials). (a) Write a computer code to evaluate the Lebesgue function L(n) (x) = n X (n) lj (x) , (3) j=0 associated to a given set of pairwise distinct nodes x0 , ..., xn . All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials) and the intellectual content of the course itself are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of Prof. Hector D. Ceniceros. 1 (b) Consider the equidistributed points xj = 1 + j(2/n) for j = 0, ..., n. Write a computer code that uses (a) to evaluate and plot L(n) (x) (evaluate L(n) (x) at a large number of points x to have a good plotting resolution, e.g. x k = 1 + k(2/ne ), k = 0, ..., ne with ne = 1000) for n = 4, 10, and 20. Estimate n for these three values of n. (c) Repeat (b) with the nodes given by xj = cos( j n ), j = 0, ..., n. Contrast the behavior of (n) L (x) and n with those corresponding to the equidistributed points in (b). 3. (a) Implement the Barycentric Formula for evaluating the interpolating polynomial for arbitrarily distributed nodes x0 , ..., xn ; you need to write a function or script that computes (n) the barycentric weights j = 1/k6=j (xj xk ) first and another code to use these values in the Barycentric Formula. Make sure to test your implementation. (b) Consider the following table of data xj 0.00 0.25 0.50 0.75 1.25 1.50 f (xj ) 0.0000 0.7071 1.0000 0.7071 -0.7071 -1.0000 Use your code in (a) to find P5 (2) as an approximation of f (2). 4. The Runge Example. Let f (x) = 1 , 1 + x2 x 2 [ 5, 5]. (4) Using your Barycentric Formula code (Prob. 3) and (5) and (6) below, evaluate and plot the interpolating polynomial of f (x) corresponding to (a) the equidistributed nodes xj = (b) (c) 5 + j(10/n), j = 0, ..., n for n = 4, 8, and 12. the nodes xj = 5 cos( j n ), j = 0, ..., n 2 Repeat (a) for f (x) = e x /5 for x 2 for n = 4, 8, 12, and 100. [ 5, 5] and comment on the result. Remark 1. It can be shown that for equidistributed nodes one can use the barycentric weights (n) j n = ( 1) , j = 0, ..., n, (5) j j where n j is the binomial coefficient (nchoosek(n,j) in Matlab). It can be shown that for the nodes xj = a+b 2 + b a 2 cos( j n ), j = 0, ..., n, in [a, b], one can use (n) j = ( 1 2( 1)j ( 1)j for j = 0 or j = n for j = 1, ..., n 1. (6) Make sure to employ (5) and (6) in your Barycentric Formula code for this problem. To plot the corresponding Pn (x) evaluate Pn (x) at a large number of points x to have a good plotting 2 resolution, e.g. x k = 5 + k(10/ne ), k = 0, ..., ne with ne = 5000. Note that your Barycentric Formula cannot be used to evaluate Pn (x) when x coincides with an interpolating node! Plot also f for comparison. Compare (a) and (b) and comment on the result in view of what you observed in Prob. 2. 3 Interpolation Hector D. Ceniceros 1 Approximation Theory Given f C[a, b], we would like to find a \"good\" approximation to it by \"simpler functions\Math 104A Homework #1 Instructor: Lihui Chai 1. Review and state the following theorems of Calculus (see for example 1.1 and 1.2 in the textbook): (a) The Intermediate Value Theorem. The intermediate value theorem states that if a continuous function, f , with an interval, [a, b], as its domain, takes values f (a) and f (b) at each end of the interval, then it also takes any value between f (a) and f (b) at some point within the interval. (b) The Mean Value Theorem. ... (c) ... ... 2. Write a computer code to implement the Composite Trapezoidal Rule quadrature \u0012 \u0013 1 1 Th [f ] = h f (x0 ) + f (x1 ) + ... + f (xN 1 ) + f (xN ) , 2 2 (1) to approximate the definite integral Z I[f ] = b f (x)dx, (2) a using the equally spaced points x0 = a, x1 = x0 + h, x2 = x0 + 2h, . . . , xN = b, where h = (b a)/N . Make sure that all your codes have a preamble which describes the purpose of the code, all the input variables, the expected output, your name, and the date of the last time you modified the code. 1 1 2 3 4 5 6 7 8 9 f u n c t i o n T = CTR( a , b ,N) % Composite T r a p e z o i d a l Rule f o r a f u n c t i o n f d e f i n e d on [ a , b ] % Input : a t h e l e f t e n d p o i n t o f t h e i n t e r v a l ; % b t h e r i g h t e n d p o i n t o f t h e i n t e r v a l ; % N t h e number o f t h e p a r t i o n used i n t h e CTR. % Output : T t h e a p p r o x i m a t i o n o f t h e i n t e g r a l g i v e n by CTR. % % Author : L i h u i Chai % Date : 2/2/17 10 11 12 13 14 15 16 17 18 19 20 h = T = T = for ( ba ) /N; 0; T + myfun ( a ) / 2 ; j = 1 : N1 x = a + j h ; T = T + myfun ( x ) ; end T = T + myfun ( b ) / 2 ; T = Th ; end % g e t t h e mesh s i z e % i n i t i a l i z e the approximation T % add t h e f i r s t node v a l u e t o T % add t h e j th node v a l u e t o T % add t h e l a s t node v a l u e t o T 21 22 23 24 25 f u n c t i o n f = myfun ( x ) % return the f u n c t i o n value at the input point x f = x exp ( x 2 ) ; end 2 3. To test your code, take f (x) = xex in [0, 1], compute the error |I[f ] Th [f ]| for h = 1/10, 1/20, 1/40, and verify that Th has a convergent trend at the expected quadratic rate. The integral I[f ] can be calculated exactly by Z I[f ] = 0 1 1 x2 1 xe dx = e = e 1. 2 0 x2 Applying the CTR with h = 1/10, 1/20, 1/40, one can obtain Th and Eh as shown in the following Table 1. As shown in the table, when h decrease, the error |Eh | decrease, and especially, as h decrease by 1/2, the error becomes to be 1/4, which indicate a quadratic convergence numerically. Table 1 h 1/10 1/20 1/40 Th |Eh | 0.8651 0.5949 102 0.8606 0.1490 102 0.8595 0.3726 103 2 x2 4. Consider the definite integral I[e 1 Z 2 ex dx. We cannot calculate its exact value but = 0 2 we can compute accurate approximations to it using Th [ex ]. Let 2 2 Th/2 [ex ] Th [ex ] q(h) = Th/4 [ex2 ] Th/2 [ex2 ] . (3) Using your code, find a value of h for which q(h) is approximately equal to 4. (a) Get an 2 2 approximation of the error, I[ex ] Th [ex ], for that particular value of h. (b) Use this error approximation to obtain the extrapolated, improved, approximation \u0011 4\u0010 2 2 2 2 Sh [ex ] = Th [ex ] + Th/2 [ex ] Th [ex ] . (4) 3 2 2 2 Explain why Sh [ex ] is more accurate and converges faster to I[ex ] than Th [ex ]. 2 Modify the CTR code to the function f (x) = ex , and compute q(h) with h = 1/10, 1/20, 1/40, 1/80, one can obtain Th and q(h) as shown in the following Table 2. As shown in the table, numerically, when h decrease, q(h) converges to 4. Especially, when h = 1/80, the difference |q(h) 4| is about 0.976 105 . Table 2 h 1/10 1/20 1/40 1/80 Th q(h) |q(h) 4| E(h) 1.46717 3.99378 0.622 103 0.452 103 1.46378 3.99844 0.156 103 0.113 103 1.46293 3.99961 0.391 104 0.283 104 1.46272 3.99990 0.976 105 0.708 105 2 2 To approximate the error Eh = I[ex ] Th [ex ], we write R(h) = 0. h0 h2 Eh = Ch2 + R(h), where lim That means R(h) is much smaller than Ch2 and the error Eh can be approximated by Ch2 . To get Ch2 , we find 2 2 I[ex ] = Th [ex ] + Ch2 + R(h) 2 = Th/2 [ex ] + C(h/2)2 + R(h/2), and then we get \u0011 4 4\u0010 2 2 Th/2 [ex ] Th [ex ] + (R(h/2) R(h)) 3 3 \u0011 4\u0010 2 2 T [ex ] Th [ex ] . 3 h/2 The approximated value of Eh calculated by the above approximation for different h can be found in Table 2. Eh Ch2 = ...... 5. ... 3 interpolate the % function f(x) = 1/(1 + x2) , x [-5,5] at equally spaced points. clear all; close all; %% f = @(x) 1./(1 + 25*x.^2); %% % Plot the original function f(x) at set of equally spaced points % figure(1); clf; x = linspace(-5,5,500); y_true = f(x); plot(x,y_true,'r','linewidth',2); hold on; %% % % Construct an interpolating of degree N at N+1 points. If we used % the barycentric form, we do not have these ill-conditioning issues! p = polyfit(xdata,ydata,N); %% % Plot the interpolating polynomial y_fit = polyval(p,x); poly_20 = plot(x,y_fit,'g','linewidth',2); plot(xdata,ydata,'k.','markersize',30); axis([-1 1 -5 5]); snapnow; %% % Let's increase one more time to see that pattern N = 40; % Increase the degree to 40 from 20. xdata = linspace(-1,1,N+1)'; ydata = f(xdata); % plot(xdata,ydata,'k.','markersize',30); % Ignore the warning about the matrix being ill-conditioned. If we used % the barycentric form, we do not have these ill-conditioning issues! p = polyfit(xdata,ydata,N); %% % Plot the interpolating polynomial y_fit = polyval(p,x); poly_40 = plot(x,y_fit,'m','linewidth',2); plot(xdata,ydata,'k.','markersize',30); axis([-1 1 -10 10]); snapnow; %% % But are we doing at least a better job in the middle? Yes! axis([-1 1 0 1]); lh = legend([poly_10, poly_20, poly_40],{'N = 10','N = 20','N = 40'}); set(lh,'fontsize',18); snapnow; % %% % % In the previous examples, our interpolating points were equally spaced % points. % % % % We can try using "Chebychev nodes". These are special nodes which allow % for a better approximation of the polynomial. % clf; y_true = f(x); plot(x,y_true,'r','linewidth',2); hold on; % Fit a 10th degree polynomial at Chebyshev nodes. N = 10; % Degree of the polynomial we try to fit t = linspace(0,pi,N+1); xdata = -cos(t); % Chebyshev nodes ydata = f(xdata); p = polyfit(xdata,ydata,N); %% % Plot the interpolating polynomial y_fit = polyval(p,x); poly_10 = plot(x,y_fit,'b','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); snapnow; % Fit a 20th degree polynomial at Chebyshev nodes. N = 20; % Degree of the polynomial we try to fit t = linspace(0,pi,N+1); xdata = -cos(t); % Chebyshev nodes ydata = f(xdata); p = polyfit(xdata,ydata,N); %% % Plot the interpolating polynomial y_fit = polyval(p,x); poly_20 = plot(x,y_fit,'g','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); snapnow; % Fit a 40th degree polynomial at Chebyshev nodes. N = 40; % Degree of the polynomial we try to fit t = linspace(0,pi,N+1); xdata = -cos(t); % Chebyshev nodes ydata = f(xdata); % Ignore the warning about the matrix being ill-conditioned. If we used % the barycentric form, we do not have these ill-conditioning issues! p = polyfit(xdata,ydata,N); %% % Plot the interpolating polynomial y_fit = polyval(p,x); poly_40 = plot(x,y_fit,'m','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); snapnow; %% % And how does the middle look? axis([-1 1 0 1]); lh = legend([poly_10, poly_20, poly_40],{'N = 10','N = 20','N = 40'}); set(lh,'fontsize',18); snapnow