Introduction to scientific Computing A Matrix Vector Approach Using Matlab Written by Charles FVan Loan 陈文斌 Wbchen(fudan. edu. cn 复日大学
Introduction to Scientific Computing -- A Matrix Vector Approach Using Matlab Written by Charles F.Van Loan 陈 文 斌 Wbchen@fudan.edu.cn 复旦大学
Chapter2 Polynomial Interpolation The vandermonde approach ° The Newton Approach Properties Special Topics
Chapter2 Polynomial Interpolation • The vandermonde Approach • The Newton Approach • Properties • Special Topics
Polynomial interpolation problem Given x…,xn( listinct)amdy…yn find a polynomialpmn(x)of degree n-1 (or less) such that pn(x)=y, for i=l: n
Polynomial interpolation problem (or less) such that p (x ) y for i :n find a polynomial p (x) of ree nGiven x ,...,x (distinct) and y ,...,y , n - i i n - n n 1 deg 1 1 1 1 1 = =
Vandermonde Approach B ase x, X pn-(x)=a+a,x+a3x+.+anx pn=1(x)=a1+02x+a3x12+…+anx 2 VI 2 2 V2 x 2 X n
Vandermonder Approach 2 1 1 1 2 3 ( ) ... − − = + + + + n n n p x a a x a x a x i n n i i i n i p x = a + a x + a x + + a x = y − − 2 1 1 1 2 3 ( ) ... = − − − − n n n n n n n n n y y y y a a a a x x x x x x x x x x x x 3 2 1 3 2 1 2 1 1 3 2 3 3 1 2 2 2 2 1 1 2 1 1 1 1 1 1 Base 1, x, x2 , … V
n= length(x);v=ones(n, n) for j=l: n n= length (x); V= ones(n, n); forⅰ=1:n Set up row i Set up row i for j=1: n V(i)=X()^j-1) V(i =(X(*ones(1, n)). c 1) end end end n=length (x) V=ones(n, n) n= length (x); V= ones(n, n)I for j=1: n for j=2: n Set up column i Set up column i for j=l: n for j=l: n V(i)=X()^j-1) V(】)=x()xV(j-1) end end end end
n = length(x); V = ones(n,n) for i=1:n % Set up row i. for j=1:n V(i,j) = x(i)^(j-1); end end n = length(x); V = ones(n,n); for i=1:n % Set up row i. V(i,:) = (x(i)*ones(1,n)).^(j- 1); end n = length(x); V = ones(n,n) for j=1:n % Set up column i. for i=1:n V(i,j) = x(i)^(j-1); end end n = length(x); V = ones(n,n) for j=2:n % Set up column i. for i=1:n V(i,j) = x(i)*V(i,j-1); end end
function a=Interp(x,y) a=Inverpv(x,y This computes the vandermonde polynomial interpolant where %o x is a column n-vector with distinct components and y is a column n-vector %o a is a column n-vector with the property that if (X)=a(1)+a(2)x+…a(n)x^(n-1) o then %p(X()=y(1),i=1n n=length(X V=ones(n, n) for i=2: n Set up column J V(:j)=X*V(:-1) end a 2n3/3 flops
function a = InterpV(x,y) % a = InverpV(x,y) % This computes the Vandermonde polynomial interpolant where % x is a column n-vector with distinct components and y is a % column n-vector. % % a is a column n-vector with the property that if % % p(x) = a(1) + a(2)x + ... a(n)x^(n-1) % then % p(x(i)) = y(i), i=1:n n = length(x); V = ones(n,n); for j=2:n % Set up column j. V(:,j) = x.*V(:,j-1); end a = V\y; 2n3 /3 flops
Nested multipication pr(=a+.+anx atx=z n-length(a); zpower=1 pval=a(1) for i=2: n zpower-z zpower pval=pval+a(i)*zpower er
Nested Multipication p x a a x x z n n = + + n = − − ( ) ... at 1 1 1 n=length(a); zpower=1; pval=a(1); for i=2:n zpower=z*zpower; pval=pval+a(i)*zpower; end
Horner's ruler m-or=a+a,x+a3x+aux =(a4x+a2)x+a2)x+a1 m=length(z) n-length(a) nF=length(a) pval=zeros(m, 1) pval=a(n) Many points Forj=1: m fori=n-1:-1:1 pval(=a(n) pval=z pval+a(1) fori=n-1:-1:1 end pval=zo *pval(+a(i) end: enc d
Horner's ruler 4 3 2 1 3 4 2 1 1 2 3 (( ) ) ( ) a x a x a x a pn x a a x a x a x = + + + − = + + + n=length(a); pval=a(n); for i=n-1:-1:1 pval=z*pval+a(i); end m=length(z);n=length(a); pval=zeros(m,1); For j=1:m pval(j)=a(n); for i=n-1:-1:1 pval(j)=z(j)*pval(j)+a(i); end;end Many points
Evalute at many points Function pval=Horner V(a, 2) mF=length(z); n=length(a) pval=a(n *ones(m, 1) fork=n-1:-1:1 pval=z "pval+a(k) en 2m flops All: 2mn flops 6 Script see next page
0 2 4 6 -2 -1 0 1 2 0 2 4 6 -2 -1 0 1 2 0 2 4 6 -2 -1 0 1 2 0 2 4 6 -2 -1 0 1 2 Evalute at many points Function pval=HornerV(a,Z) m=length(z);n=length(a); pval=a(n)*ones(m,1); for k=n-1:-1:1 pval=z.*pval+a(k); end 2m flops All:2mn flops Script see next page
o Script File: Show V 1% Plots 4 random cubic interpolants of sin(x)on [0, 2pi] %o Uses the vandermonde method close a :0-linspace(0, 2 pi, 100) yo= sin(x0) for eg=1: 4 Display cubic interpolants of sin(x). The abscissas are chosen x=2*pi*sort(rand( 4, 1) randomly sin(X a=Interp( y) p Val= Horner v(a,]O) subplot(2, 2, eg plot( xo, y0, 0, p Val, x,y, * axiS([02*p-22]) end
% Script File: ShowV % Plots 4 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:4 x = 2*pi*sort(rand(4,1)); y = sin(x); a = InterpV(x,y); pVal = HornerV(a,x0); subplot(2,2,eg) plot(x0,y0,x0,pVal,'--',x,y,'*') axis([0 2*pi -2 2]) end Display cubic interpolants of sin(x). The abscissas are chosen randomly
