Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 6 Prof. Wickerhauser Exercise 2(b) of Section 5.1, p.259. MATLAB CODE: <>x=[-6 -2 0 2 6]; y=[-5.3 -3.5 -1.7 0.2 4]; f=[-6 -2.84 -1.26 0.32 3.48]; <>one=[1 1 1 1 1]; n=5; <>xmean=one*x'/n XMEA = 0. <>ymean=one*y'/n YMEA = -1.2600 <>x-xmean*one ANS = -6. -2. 0. 2. 6. <>sxy = (x-xmean*one)*(y-ymean*one)' SXY = 63.2000 <>sx2 = (x-xmean*one)*(x-xmean*one)' SX2 = 80. <>a = sxy/sx2 A = .7900 <>b=ymean-a*xmean B = -1.2600 <>l=a*x+b*one L = -6.0000 -2.8400 -1.2600 .3200 3.4800 <>l-f ANS = 1.0E-15 * .0000 .0000 -.2220 -.3886 -.8882 <>l-y ANS = -.7000 .6600 .4400 .1200 -.5200 <>e2=(l-y)*(l-y)' E2 = .5299 Exercise 5 of Section 5.1, p.260. Let SX2, SX1, SY1, and SXY be the sums of x^2, x, y, and x*y, respectively, from Equation Eq. 10 of p.255. Then Eq. 10 can be written in matrix form as: / SX2 SX1 \ /A\ = /SXY\ \ SX1 N / \B/ \SY1/ This 2x2 system can be solved explicitly: /A\ = -1 / N -SX1 \ /SXY\ , \B/ D \ -SX1 SX2 / \SY1/ where D = SX2*N-SX1*SX1. Multiplying this out gives A = ( N*SXY - SX1*SY1)/D; B = (SX2*SY1 - SX1*SXY)/D Exercise 6 of Section 5.1, p.260. Let xm = SX1/N be the mean value of x(1),...,x(N). Then note that [x(k)-xm]*[x(k)-xm] = x(k)*x(k) - 2*xm*x(k) + xm*xm, so summing this over k=1,2,...,N gives Sum [x(k)-xm]*[x(k)-xm] = SX2 - 2*xm*SX1 + N*xm*xm k = SX2 - 2*SX1*SX1/N + SX1*SX1/N = SX2 - SX1*SX1/N = D/N This shows that D = N*Sum [x(k)-xm]*[x(k)-xm], k which is nonzero if x(1),...,x(N) contains at least 2 distinct points, since then [x(k)-xm] must be nonzero for some k. Algorithm 4(a,b,c) of Section 5.1, p.262. MATLAB CODE: OUTPUT: % generate the function: n=50; x=[1:n] * 0.1; y=x+cos(sqrt([1:n])); %(a) compute least-squares line y = ax+b: sx1=sum(x); sy1=sum(y); sx2=sum(x.*x); sxy=sum(x.*y); d=sx2*n-sx1*sx1; a=( n*sxy - sx1*sy1)/d A = 1.4186 b=(sx2*sy1 - sx1*sxy)/d B = -.8808 l=a*x+b*ones(x); %(b) compute the RMS error: e2=norm(l-y)/sqrt(n) E2 = .4032 %(c) plot the error: plot(l-y) * ** ** * * * * * * * * * * * * * * ** * * * * * ** ** ** * ** * * * ** ** * ** * * * * Exercise 2(c) of Section 5.2, p.275. MATLAB CODE: OUTPUT: x=[-2 -1 0 1 2]; y=[10 1 0 2 9]; n=5; % compute least-squares parabols y = ax^2+bx+c: sx1=sum(x); sx2=sum(x.*x); sx3=sum(x.*x.*x); sx4=sum(x.*x.*x.*x); sy1=sum(y); sxy=sum(x.*y); sxxy=sum(x.*x.*y); mx=[sx4 sx3 sx2 ; sx3 sx2 sx1 ; sx2 sx1 n]; rh=[sxxy sxy sy1]' ; coef=mx\rh; a=coef(1), A = 2.5000 b=coef(2), B = -.1000 c=coef(3), C = -.6000 p=a*x.*x+b*x+c*ones(x), P = 9.6000 2.0000 -.6000 1.8000 9.2000 y, Y = 10. 1. 0. 2. 9. y-p, = .4000 -1.0000 .6000 .2000 -.2000 Another way is to use the degree-2 general expression for least-squares polynomials: <> x=[-2 -1 0 1 2]' ; y=[10 1 0 2 9]' ; f=[ones(x) x x.*x] ; <> c=(f'*f)\(f'*y) C = -.6000 -.1000 2.5000 <> (f*c)' ANS = 9.6000 2.0000 -.6000 1.8000 9.2000 <> e2=norm(f*c-y)/sqrt(5) E2 = .5657 Exercise 12 of Section 5.2, p.277. Linearization is achieved by the change of variable Y=y, X=ln(x). The new variables X,Y are related by Y=AX+B, if the old variables are related by y=A ln(x) + B. Exercise 2(a,b) of Section 5.3, p.294. (a) If S(x)=a0+a1*x +a2*x^2 +a3*x^3, then S'(x)=a1+2*a2*x +3*a3*x^2. Substituting x=1 and x=2 and the given values gives the stated relations. (b)MATLAB CODE: OUTPUT: s=[ 1 1 1 1 ; 0 1 2 3 ; 1 2 4 8 ; 0 1 4 12 ]; rh=[3 -4 1 2]' ; a=s\rh, A = 5. 2. -6. 2. x=[0:30]*0.1; p=a(1)*ones(x)+x.*(a(2)*ones(x)+x.*(a(3)*ones(x)+x*a(4))); plot(p) * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Algorithm 5 of Section 5.3, p.297. First modify the clamped cubic spline program to produce the natural cubic spline: function S=nsfit(X,Y) %Input - X is the 1xn abscissa vector % - Y is the 1xn ordinate vector %Output - S: rows of S are the coefficients for the cubic interpolants % NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248-2 % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458 N=length(X)-1; H=diff(X); D=diff(Y)./H; A=H(2:N-1); B=2*(H(1:N-1)+H(2:N)); C=H(2:N); C=H(2:N); U=6*diff(D); %Natural spline endpoint constraints % B(1)=B(1)-H(1)/2; % no adjustment needed % U(1)=U(1)-3*(D(1)-dx0); % no adjustment needed % B(N-1)=B(N-1)-H(N)/2; % no adjustment needed % U(N-1)=U(N-1)-3*(dxn-D(N)); % no adjustment needed for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end % Natural spline endpoint constraints M(1)=0; M(N+1)=0; for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end Next, modify the code on page 294 to generate points to plot on the piecewise cubic spline: % Evaluation of cubic splines on M subintervals % [X(1),X(2)], [X(2),X(3)],..., [X(M),X(M+1)] function [Xout,Yout]=seval(S,X,steps) M=length(X)-1; Xout = X(1); % initialize output abscissas Yout = S(1,4); % initialize output ordinates for k=1:M delX = (X(k+1)-X(k))/steps; % increment per step, kth subinterval newX = X(k):delX:X(k+1); % abscissas, kth subinterval newY = polyval(S(k,:), newX-X(k)); % ordinates, kth subinterval Xout = [Xout, newX]; Yout = [Yout, newY]; end Finally, compute the integral on [xk,x{k+1}] of the cubic spline S(x) = s0k + s1k(x-xk) + s2k(x-xk)^2 + s3k(x-xk)^3 by elementary methods: Int = s0k*dx +(1/2)*s1k*dx^2 + (1/3)*s1k*dx^3 + (1/4)*s1k*dx^4 , where dx = x{k+1}-xk. Since dx=1 for all k=1,...,12 in this problem, the integral of the cubic spline is given by Int = s0k +(1/2)*s1k + (1/3)*s1k + (1/4)*s1k However, the function nsfit() returns the spline coefficients in reverse order: s0k <----> S(k,4) s1k <----> S(k,3) s2k <----> S(k,2) s3k <----> S(k,1) so the Matlab formula for the integral is >> sum(S*[1/4 1/3 1/2 1]') ans = 656.9545 and the average is computed by the following Matlab code: >> parts = S*[1/4 1/3 1/2 1]'; avg = sum(parts)/length(parts) avg = 59.7231 Compare this with the average of the hourly temperatures: >> sum(Y)/length(Y) ans = 60 INPUT/OUTPUT: >> X=1:12; Y=[58 58 58 58 57 57 57 58 60 64 67 68]; >> plot(X,Y) % piecewise linear plot of the data >> S=nsfit(X,Y) S = -0.0247 0 0.0247 58.0000 0.1236 -0.0742 -0.0494 58.0000 -0.4697 0.2966 0.1730 58.0000 0.7551 -1.1124 -0.6427 58.0000 -0.5509 1.1530 -0.6021 57.0000 0.4483 -0.4996 0.0513 57.0000 -0.2425 0.8454 0.3971 57.0000 0.5216 0.1180 1.3605 58.0000 -0.8438 1.6827 3.1611 60.0000 -0.1464 -0.8487 3.9951 64.0000 0.4293 -1.2878 1.8585 67.0000 >> [Xout,Yout]=seval(S,X,50); >> plot(Xout,Yout) % cubic spline fit to the data >> plot(Xout,Yout,X,Y) % spline with original data See the plot at ./5-3-a5.pdf Exercise 2 of Section 5.4, p.307. The function f is even --- f(-x)=f(x) --- so its Fourier series contains only cosine functions: f(x) = a(0)/2 + Sum(n=1,2,...) a(n) cos(nx), where a(n)=(1/pi)Integral(-pi,pi) f(x) cos(nx) dx for n=0,1,2,... These integrals may be computed exactly, using the calculus: pi*a(n) = Integral(-pi,0)(pi/2 + x) cos(nx) dx + Integral(0,pi)(pi/2 - x) cos(nx) dx = Integral(0,pi)(pi/2 - x) cos(nx) dx + Integral(0,pi)(pi/2 - x) cos(nx) dx = 2 Integral(0,pi)(pi/2 - x) cos(nx) dx, where we have made the substitution x <-- -x in the first integral and then reduced and combined the results. Now, if n=0, we get pi*a(0) = 2 Integral(0,pi)(pi/2 - x) dx = 0. Otherwise, for n>0, we have can evaluate the integral in two pieces involving cosine. The first piece is trivial: 2 Integral(0,pi) (pi/2) cos(nx) dx = 2 pi [sin(n pi)/n - sin(0)/n] = 0; The second can be integrated by parts, using u=x and dv=cos(nx)dx. This yields du=dx and v=sin(nx)/n, giving: -2 Integral(0,pi) x cos(nx) dx = [x sin(n x)/n]|(0,pi) - Integral(0,pi) sin(nx)/n dx = (-2) { [0 - 0] + cos(n pi)/n^2 - 1/n^2 }, which is 0 if n is even and 4/n^2 if n is odd. Thus a(n)=0 for even n and a(n)=4/[pi*n^2] if n is odd, so f(x) = Sum(n=1,3,5,...) 4*cos(n x)/(pi*n^2) Exercise 7 of Section 5.4, p.307. If we set x=0 in Exercise 2, the function f(x) evaluates to pi/2, while its Fourier series simplifies to Sum(n=1,3,5,...) 4/(pi n^2). Combining constants yields the identity pi^2 / 8 = Sum(n=1,3,5,...) 1/n^2. Algorithm 1 of Section 5.5, p.319. The following commands implement the generation of samples along a Bezier curve of degree N from an (N+1)x2 input matrix of N+1 control points: function [x,y]=bezier(M,dt) %Input - M is an (N+1)x2 matrix of N+1 control points % - dt is the increment in the parameterization by 0> Ma = [1 3; 3 -1; 2 4; 3 0]; [x y]=bezier(Ma,0.01); plot(x,y) >> Mb = [-2 3; -1 3; 3 5; 3 4; 2 3]; [x y]=bezier(Mb,0.01); plot(x,y) >> Mc=[1 1;2 2;3 4;4 4;5 2;6 1]; [x y]=bezier(Mc,0.01); plot(x,y) The results are viewable in files ./5-5-a2a.pdf, ./5-5-a2b.pdf, and ./5-5-a2c.pdf For the next two problems, set cas(x)= cos(2*PI* x)+sin(2*PI*x). Problem H06.1. Fix N and show that, if n and m are in the range [0,1,...,N-1] with n not equal to m, then the sum from k=0 to k=N-1 of cas(n*k/N)*cas(m*k/N) equals 0. What is the value when n=m? Fix N. Then for all n,m, using trigonometric identities we show that Sum(k=0,...,N-1) cas(n*k/N) cas(m*k/N) = Sum(k=0,...,N-1) [ cos(2 pi n k / N) cos(2 pi m k / N) + cos(2 pi n k / N) sin(2 pi m k / N) + sin(2 pi n k / N) cos(2 pi m k / N) + sin(2 pi n k / N) sin(2 pi m k / N) ] = Sum(k=0,...,N-1) [ cos(2 pi(n-m)k/N) + sin(2 pi(n+m)k/N)]. We can analyze the two summands individually. First note that cos x = Re exp(ix). Thus, for n != m, we have exp(2i pi(n-m)k/N) != 1, and we can use the geometric sum formula: Sum(k=0,...,N-1) cos(2 pi(n-m)k/N) = Re Sum(k=0,...,N-1) exp(2i pi(n-m)k/N) = Re [(1-exp(2i pi(n-m)k)]/[1-exp(2i pi(n-m)k/N)] = 0. Second, note that sin x = Im exp(ix). If n+m=0, then the second summand is obviously 0 since all its terms are zero. Otherwise, we again apply the geometric sum formula to get Sum(k=0,...,N-1) sin(2 pi(n+m)k/N) = Im Sum(k=0,...,N-1) exp(2i pi(n+m)k/N) = Im [1-exp(2i pi(n+m)k)]/[1-exp(2i pi(n+m)k/N)] = 0. This shows that the inner product of two ``cas'' functions is zero if n != m. If n=m, then the second summand is still 0 but the first becomes 1 for all k and thus the inner product is N. Problem H06.2. Show that, for N=2M and any function f=f(k) defined on 0,1,...,N-1, we have that the sum from k=0 to k=N-1 of f(k)*cas(n*k/N) equals the sum from k=0 to k=M-1 of f(2*k)*cas(n*k/M) plus the sum from k=0 to k=M-1 of f(2*k+1)*cas(n*k/M + n/2M). For N=2M and any function f=f(k) defined on 0,1,...,N-1, Sum(m=0,...,2M-1) f(m) cas(m k/2M) = Sum(m=0,...,M-1) f(2m) cas(n(2m)/2M) +Sum(m=0,...,M-1) f(2m+1) cas(n(2m+1)/2M) = Sum(m=0,...,M-1) f(2m) cas(nm/M) +Sum(m=0,...,M-1} f(2m+1) cas(nm/M+ n/2M).