Condsplinemat18 ans 32707 recall the vandermonde

This preview shows page 39 - 42 out of 60 pages.

>cond(splinemat(18)) ans = 32.707 Recall the Vandermonde matrix had a condition number of 1.8822e+14 . This shows that the system of equations for the splines is very much better conditioned, by 13 orders of magnitude!! Code for splinemat.m and plotspline.m function S=splinemat(n) L=[1 1 1;3 2 1;6 2 0]; M=[0 0 0;0 0 -1; 0 -2 0]; Z=zeros(3,3); T=[0 0 0;0 2 0; 0 0 0]; V=[1 1 1;0 0 0;6 2 0]; S=zeros(3*(n-1),3*(n-1)); for k=[1:n-2] for l=[1:k-1] 39
S(3*k-2:3*k,3*l-2:3*l) = Z; end S(3*k-2:3*k,3*k-2:3*k) = L; S(3*k-2:3*k,3*k+1:3*k+3) = M; for l=[k+2:n-1] S(3*k-2:3*k,3*l-2:3*l) = Z; end end S(3*(n-1)-2:3*(n-1),1:3)=T; for l=[2:n-2] S(3*(n-1)-2:3*(n-1),3*l-2:3*l) = Z; end S(3*(n-1)-2:3*(n-1),3*(n-1)-2:3*(n-1))=V; end function plotspline(X,Y) n=length(X); L=X(2)-X(1); S=splinemat(n); b=zeros(1,3*(n-1)); for k=[1:n-1] b(3*k-2)=Y(k+1)-Y(k); b(3*k-1)=0; b(3*k)=0; end a=S\b’; npoints=50; XL=[]; YL=[]; for k=[1:n-1] XL = [XL linspace(X(k),X(k+1),npoints)]; p = [a(3*k-2),a(3*k-1),a(3*k),Y(k)]; XLL = (linspace(X(k),X(k+1),npoints) - X(k)*ones(1,npoints))/L; YL = [YL polyval(p,XLL)]; end plot(X,Y,’o’) hold on 40
plot(XL,YL) hold off I.2.6. The linear equations for cubic splines: version 2 This version follows Numerical Recipes in C, by Press, et. al. Given points ( x 1 , y 1 ) , . . . , ( x n , y n ) with x 1 < x 2 < · · · < x n we wish to find a collection of cubic polynomials p 1 ( x ) , p 2 ( x ) , . . . , p n - 1 ( x ) such that p i ( x ) is our interpolating function in the interval [ x i , x i +1 ]. We require that 1. p i ( x i ) = y i and p i ( x i +1 ) = y i +1 for i = 1 , . . . , n - 1 (These are 2( n - 1) equations stating that each p i has the correct values at the left and right endpoints of its interval.) 2. p 0 i ( x i +1 ) = p 0 i +1 ( x i +1 ) for i = 1 , . . . , n - 2 (These are n - 1 equations stating that first derivatives are continuous.) 3. p 00 i ( x i +1 ) = p 00 i +1 ( x i +1 ) for i = 1 , . . . , n - 2 (These are n - 1 equations stating that second derivatives are continuous.) 4. p 00 1 ( x 1 ) = p 00 n - 1 ( x n ) = 0 (These are 2 equations stating that the endpoint conditions hold.) There are a total of 2( n - 1) + 2( n - 2) + 2 = 4 n - 4 equations. Each polynomial p i can be written in the form A i x 3 + B i x 2 + C i x + D i and so contains 4 unknown coefficients. Since there are n - 1 polynomials, this gives a total of 4( n - 1) = 4 n - 4 unknowns. The equations above form a system 4 n - 4 linear equations in 4 n - 4 unknown coefficients. These have a unique solution which can be found and plotted using MATLAB/Octave. Even though it turns out to be a bad idea to carry this out directly, we will go ahead anyway so we can compare this version with the more efficient one coming up next. Let’s write down the equations for the case of n = 3 points ( x 1 , y 1 ) , ( x 2 , y 2 ) , ( x 3 , y 3 ). In this case we want to find two polynomials, p 1 ( x ) = A 1 x 3 + B 1 x 2 + C 1 x + D 1 and p 2 ( x ) = A 2 x 3 + B 2 x 2 + C 2 x + D 2 so there are eight unknowns A 1 , B 1 , C 1 , D 1 , A 2 , B 2 , C 2 , D 2 . The equations are 1 . p 1 ( x 1 ) = y 1 A 1 x 3 1 + B 1 x 2 1 + C 1 x 1 + D 1 = y 1 p 1 ( x 2 ) = y 2 A 1 x 3 2 + B 1 x 2 2 + C 1 x 2 + D 1 = y 2 p 2 ( x 2 ) = y 2 A 2 x 3 2 + B 2 x 2 2 + C 2 x 2 + D 2 = y 2 p 2 ( x 3 ) = y 3 A 2 x 3 3 + B 2 x 2 3 + C 2 x 3 + D 2 = y 3 2 . p 0 1 ( x 2 ) = p 0 2 ( x 2 ) 3 A 1 x 2 2 + 2 B 1 x 2 + C 1 - 3 A 2 x 2 2 - 2 B 2 x 2 - C 2 = 0 3 . p 00 1 ( x 2 ) = p 00 2 ( x 2 ) 6 A 1 x 2 + 2 B 1 - 6 A 2 x 2 - 2 B 2 = 0 4 . p 00 1 ( x 1 ) = 0 6 A 1 x 1 + 2 B 1 = 0 p 00 2 ( x 3 ) = 0

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture