lagrange插值分段线性插值matlab代码

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

输出结果：

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)

输出结果：

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个点，其中，54个有效

首先保存你一个M文件：

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

axcs('position,,[0011])

[x,y]=ginput;

然后在commandwindow里,输入以下内容：

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,'，)；

对应的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个点，全部有效

首先保存一个M文件：

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

axcs('position,,[0011])

[x,y]=ginput;

然后在commandwindow里,输入以下内容：

n=length(x);

s=(l:n),;

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

u=pchip(s,x,t)；

v=pchip(s,y,t);

elfreset

对应的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的，这个能将之间画出弧度，而pchip更像是直接用线段将点依

次连接得到的。

使用的是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?

解：

⑶

首先保存一个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里输入以卜内容:

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时，图像：

n=10时，图像:

n=100时，图像:

n=1000时，图像:

可以看出，将卜1,1］做n・l等分的n个插值点，在卜0.92,1)的区间内，随着n趋近于8时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

输出结果：

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时,

0.4121-0.90770.3897

-0.9363-0.8057-0.3658

-0.44620.64690.9004

n=10时，

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时，数据太大，没运行出来。

可以看出，将卜1,1］做n-1等分的n个插值点，在卜0.92,0.92］的区间内，随着n趋近于8时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)

输出结果：

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