Lagrange插值：

x=0:3;

y=[-5,-6z-l,16];

n=length(x);

symsq;

fork=l:n

fenmu=l;

p=l；

forj=11n

if(j-=k)

fenmu=fennu*(x(k)-x(j))

p=conv(pzpoly(x(j)))

end

end

c(k,:)=p*y(k)/fenmu

end

a=zeros(1,n);

fori-l:n

forj=1:n

a(i)=a(i)+c(j,i)

end

end

輸出結(jié)果：

fenmu=

-1

P=

1-1

fenmu=

2

P=

1-32

fenmu=

-6

P=

1-611-6

c=

0.8333-5.00009.1667-5.0000

fenmu=

1

P=

10

fenmu=

-1

p=

1-20

fenmu=

2

P=

1-560

c=

0.8333-5.00009.1667-5.0000

-3.000015.0000-18.00000

fenmu=

2

P=

10

fenmu=

2

P=

1-10

fenmu=

-2

P=

1-430

c=

0.8333-5.00009.1667-5.0000

-3.000015.0000-18.00000

0.5000-2.00001.50000

fenmu=

3

P=

10

fenmu=

6

P=

1-10

fenmu=

6

P=

1-320

c=

0.8333-5.00009.1667-5.0000

-3.000015.0000-18.00000

0.5000-2.00001.50000

2.6667-8.00005.33330

a=

0.8333000

a=

-2.1667000

a=

-1.6667000

a=

1000

a=

1-500

a=

11000

3=

1800

a=

1000

a=

1.000009.16670

a=

1.00000-8.83330

a=

1.00000-7.33330

a=

1.00000-2.00000

a=

1.00000-2.0000-5.0000

a=

1.00000-2.0000-5.0000

a=

1.00000-2.0000-5.0000

a=

1.00000-2.0000-5.0000

分段線性插值：

先保存M文件：

x=l:6;

y=[7168251224],

u=5.3;

delta=diff(y)./diff(x);

n=length(x);

forj=2:(n-1)

ifx(j)<u

k=j;

end

end

在commandwindow中輸入:

s=u-x(k);

v=y(k.)+s.*delta(k)

輸出結(jié)果：

v=

15.6000

3.4.Makeaplotofyoiirhand.Startwith

figure(position),get(0,'screensize'))

axes('position',[0Oil])

[x,y]=ginput;

Placeyourhandonthecomputerscreen.Usethemousetoselectafew

dozenpointsoutliningyourhand.Terminatetheginpurwithacarriage

return.Youmightfinditeasiertotraceyourhandonapieceofpaperand

thenputthepaperonthecomputerscreen.Youshouldl>eabletoseethe

ginputcursorthroughthepaper.(Savethesedata.Wewillrefertothemin

otherexerciseslaterinthisbook.)

Figure3.11.Ahand.

Nowthinkofxandyastwofunctionsofanindependentvariablethatgoes

fromonetothenumberofpointsjroucollected.Youcaninterpolateboth

functionsonafinergridandplottheresultwith

n=lengrh(x);

s=(1:n),;

z?(1:.05:n)1;

u=splinetx(s,x,t);

v=splinetx(s,y,t);

elfreset

plot(x,y,1',u,v,;

Dothesamethingwithpchiptx.Whichdoyouprefer?

Figure3.11istheplotofmyhand.Canyoutellifitwasdonewithsplinetx

orpchiptx?

解：

第一種做法，用spline,共55個(gè)點(diǎn)，其中，54個(gè)有效

首先保存你一個(gè)M文件：

figure('position',get(0,'screensize'))

axcs('position,,[0011])

[x,y]=ginput;

然后在commandwindow里,輸入以下內(nèi)容：

n=length(x);

s=(l:n),;

t=(l:.O5:n),;

u=spline(s,x,t);

v=splinc(s,y,t);

elfreset

plot(x,y,'.',u,v,'，)；

對(duì)應(yīng)的x、y值:

0.35729170.2536145

0.35729170.2909639

0.35034720.3403614

0.34618060.4259036

0.34270830.5271084

0.32534720.6162651

0.30659720.6873494

0.2906250.7524096

0.28923610.7933735

0.29548610.796988

0.32256940.7548193

0.3406250.6849398

0.36909720.6150602

0.38645830.6126506

0.38993060.7259036

0.39270830.8066265

0.39201390.8993976

0.40243060.9295181

0.42395830.8933735

0.42395830.8078313

0.42951390.7343373

0.43159720.6451807

0.44409720.6439759

0.45659720.7439759

0.47048610.8451807

0.47673610.9054217

0.49618060.9463855

0.50868060.876506

0.50451390.818G747

0.50104170.7524096

0.48923610.6403614

0.5031250.6295181

0.50520830.6271084

0.53229170.7090361

0.55104170.763253

0.57395830.8355422

0.59618060.8572289

0.59479170.7837349

0.57534720.7090361

0.55798610.6391566

0.53576390.5668675

0.53229170.5283133

0.53506940.4789157

0.5656250.536747

0.59479170.5933735

0.62534720.610241

0.63229170.5728916

0.6156250.5331325

0.60034720.4993976

0.57881940.4415663

0.5593750.3716867

0.52951390.2957831

0.49756940.2403614

0.47118060.2018072

0.66076390.3090361

第二種做法，用pchip,共52個(gè)點(diǎn)，全部有效

首先保存一個(gè)M文件：

figure('position',get(0,'screensize'))

axcs('position,,[0011])

[x,y]=ginput;

然后在commandwindow里,輸入以下內(nèi)容：

n=length(x);

s=(l:n),;

t=(l:.O5:n),;

u=pchip(s,x,t)；

v=pchip(s,y,t);

elfreset

對(duì)應(yīng)的x、y值:

0.51909720.8487952

0.50520830.7512048

0.49479170.6789157

0.51006940.6692771

0.53993060.7355422

0.57534720.8174699

0.5968750.8620482

0.61909720.8777108

0.61493060.8138554

0.58784720.7427711

0.58784720.7427711

0.56354170.6716867

0.53506940.603012

0.5281250.563253

0.5281250.5259036

0.5656250.5801205

0.60520830.6271084

0.6343750.6186747

0.61909720.5716867

0.58784720.523494

0.53645830.4126506

0.49618060.3210843

0.4593750.2753012

我更喜歡第一種，用spline的，這個(gè)能將之間畫出弧度，而pchip更像是直接用線段將點(diǎn)依

次連接得到的。

使用的是splinetx。

3.9.TheM-filerungeinterp.mprovidesanexperimentwithafamouspolynomial

interpolationproblemduetoCarlRunge.Let

f3=TT—

an<lletPn(x)denotethepolynomialofdegreen—1thatinterpolatesf(x)at

nequallyspacedpointsontheinterval-1<i<1.Rmigeaskedwhether,

asnincreases.Pn(x)convergesto/(□").Theanswerisyesforsomei,but

noforothers.

(a)ForwhatxdocsPn(x)—*/(z)asn—?oo?

(b)Changethedistributionoftheinterpolationpointssothattheyarenot

equallyspaced.HDWdocsthisaffectconvergence?Canyoufindadistribution

sothatPn(x)—*f(x)foralla*intheinterval?

解：

⑶

首先保存一個(gè)M文件：

n=3;

xishu=2/(n-1);

x=-l:xishu:1;

y=l./(1+25.*x.*x);

fork=l:n

fenmu=l;

p=l；

forj=l:n

if(j~=k)

fenmu=fennu*(x(k)-x(j));

p=conv(p,poly(x(j)));

end

end

c(k,:)=p*y(k)/fenmu;

end

a=zeros(1,n);

fori=l:n

forj=l:n

a(i)=a(i)+c(j,i);

end

pnd

然后在commandwindow里輸入以卜內(nèi)容:

plot(x,y；*')；

holdon;

plotfx^；*');

holdon;

xl=-l:0.01:l;

yi=o；

fori=l:n

yl=yl.*xl+a(i)

end

y2=l./(l+25.*xl.*xl);

plot(xl,yl/b')；

holdon;

plot(xl,y2；g');

即有n=3時(shí)，圖像：

n=10時(shí)，圖像:

n=100時(shí)，圖像:

n=1000時(shí)，圖像:

可以看出，將卜1,1］做n?l等分的n個(gè)插值點(diǎn)，在卜0.92,1)的區(qū)間內(nèi)，隨著n趨近于8時(shí)P￡x)

趨近于f(x)o

(b)

先保存M文件：

n=2;

x=2.*rand(n)-1

y=l./(1+25.*x.*x);

n=n八2;

%lagrangc?????i§

fork=l:n

fenmu=l;

p=l;

forj=1:n

if(j~=k)

fADTnu=fpnnn*(x(k)-x(j));

p=conv(p,poly(x(j)));

end

end

c(k,:)=p*y(k)/fenmu;

end

a=zeros(1,n);

fori=l:n

forj=l:n

a(i)-a(i)Ic(j,i);

end

end

輸出結(jié)果：

x=

0.9150-0.6848

0.92980.9412

然后在commandwindow里輸入：

plot(x,y；*');

holdon;

xl=-l:0.01:l;

yi=o；

fori=l:n

yl=yl.*xl+a(i)

end

y2=l./(l+25.*xl.*xl);

plot(xl,yl；b');

holdon;

plot(xlzy2/g')；

得到以下幾幅圖：

igure1-X

Eil?Edit乂ievInsert工ooIsfiesktop{indovHdLp

"U)Q?、-、門S4.乂?Q□Id-□

1??<<IIII

*

0.075--

0.07--

0.066--

0.06--

0.055--

0.05--

0.045■%-

*

0.041----------1----------1----------1----------1----------1----------1----------1----------1----------

飛.8-0.6-04-020020.40.60.81

n=3時(shí),

0.4121-0.90770.3897

-0.9363-0.8057-0.3658

-0.44620.64690.9004

n=10時(shí)，

x=

Columns1through7

-0.32100.9661-0.65780.71100.1660-0.7845-0.6425

0.9033-0.3971-0.93480.2895-0.49640.8126-0.1542

0.84070.40220.1224-0.2475-0.41910.7593-0.8115

-0.89460.33270.7637-0.61820.23420.63550.1970

0.47570.07830.3384-0.1435-0.4694-0.4785-0.0582

-0.46180.3962-0.6191-0.03600.64880.18870.3919

-0.15430.3331-0.2622-0.75880.9653-0.95500.3998

0.0957-0.6437-0.07850.17900.4605-0.14950.2771

0.8855-0.74400.9633-0.5476-0.3122-0.3746-0.9328

-0.16450.9982-0.6872-0.23080.1681-0.6770-0.8624

Columns8through10

-0.36080.22190.7507

0.06170.55760.0361

0.3089-0.15310.8872

-0.1848-0.81840.2754

0.6400-0.46710.9154

0.4367-0.6927-0.5186

0.9373-0.43800.3522

0.0627-0.1198-0.4219

-0.34970.05430.3436

-0.7887-0.08520.3903

n=10時(shí)，數(shù)據(jù)太大，沒運(yùn)行出來。

可以看出，將卜1,1］做n-1等分的n個(gè)插值點(diǎn)，在卜0.92,0.92］的區(qū)間內(nèi)，隨著n趨近于8時(shí)P￡x)

趨近于f(X)o

3.18.(a)Ifyouwanttointerpolatecensusdataontheinterval1900<t<2000

withapolynomial,

109

P(t)=cjt+c2rH-----bc10t4-cH,

youmightbetemptedtousetheVandermondematrixgeneratedby

t=1900:10:2000

V=vander(t)

Whyisthisareallybadidea?

(b)Investigatecenteringandscalingtheindependentvariable.Plotsome

data,pulldowntheToolsmenuonthefigurewindow,selectBasicFitting,

andfindthecheckboxaboutcenteringandscaling.Whatdoesthischeck

lx)xdo?

(c)Replacethevariabletwith

t-

s=-----”

a

ThisleadstoamodifiedpolynomialP(s).Howareitscoefficientsrelatedto

thoseof尸(￡)?WhathappenstotheVandermondematrix?Whatvaluesof

〃andaleadtoareasonablywellcoiuiitionedVandernioiidematrix?One

possibilityis

mu-mean(t)

sigma=srd(r)

butaretherebettervalues?

解：

(a)

t=1900:10:2000

V=vander(t)

輸出結(jié)果：

t=

Columns1through6

190019101920193019401950

Columns7through11

19601970198019902000

V=

1.0e+033*

Columns1through7

0.61310.00030.00000.00000.00000.00000.0000

0.64620.00030.00000.00000.00000.00000.0000

0.68080.000