lagrange插值分段线性插值matlab代码.docx
《lagrange插值分段线性插值matlab代码.docx》由会员分享,可在线阅读,更多相关《lagrange插值分段线性插值matlab代码.docx(25页珍藏版)》请在冰点文库上搜索。
lagrange插值分段线性插值matlab代码
Lagrange插值:
x=0:
3;
y=[-5,-6,-1,16];
n=length(x);
symsq;
fork=1:
n
fenmu=1;
p=1;
forj=1:
n
if(j~=k)
fenmu=fenmu*(x(k)-x(j))
p=conv(p,poly(x(j)))
end
end
c(k,:
)=p*y(k)/fenmu
end
a=zeros(1,n);
fori=1:
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
a=
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=1:
6;
y=[7168251224];
u=5.3;
delta=diff(y)./diff(x);
n=length(x);
forj=2:
(n-1)
ifx(j)
k=j;
end
end
在commandwindow中输入:
s=u-x(k);
v=y(k)+s.*delta(k)
输出结果:
v=
15.6000
解:
第一种做法,用spline,共55个点,其中,54个有效
首先保存你一个M文件:
figure('position',get(0,'screensize'))
axes('position',[0011])
[x,y]=ginput;
然后在commandwindow里输入以下内容:
n=length(x);
s=(1:
n)';
t=(1:
.05:
n)';
u=spline(s,x,t);
v=spline(s,y,t);
clfreset
plot(x,y,'.',u,v,'-');
对应的x、y值:
0.3572917
0.2536145
0.3572917
0.2909639
0.3503472
0.3403614
0.3461806
0.4259036
0.3427083
0.5271084
0.3253472
0.6162651
0.3065972
0.6873494
0.290625
0.7524096
0.2892361
0.7933735
0.2954861
0.796988
0.3225694
0.7548193
0.340625
0.6849398
0.3690972
0.6150602
0.3864583
0.6126506
0.3899306
0.7259036
0.3927083
0.8066265
0.3920139
0.8993976
0.4024306
0.9295181
0.4239583
0.8933735
0.4239583
0.8078313
0.4295139
0.7343373
0.4315972
0.6451807
0.4440972
0.6439759
0.4565972
0.7439759
0.4704861
0.8451807
0.4767361
0.9054217
0.4961806
0.9463855
0.5086806
0.876506
0.5045139
0.8186747
0.5010417
0.7524096
0.4892361
0.6403614
0.503125
0.6295181
0.5052083
0.6271084
0.5322917
0.7090361
0.5510417
0.763253
0.5739583
0.8355422
0.5961806
0.8572289
0.5947917
0.7837349
0.5753472
0.7090361
0.5579861
0.6391566
0.5357639
0.5668675
0.5322917
0.5283133
0.5350694
0.4789157
0.565625
0.536747
0.5947917
0.5933735
0.6253472
0.610241
0.6322917
0.5728916
0.615625
0.5331325
0.6003472
0.4993976
0.5788194
0.4415663
0.559375
0.3716867
0.5295139
0.2957831
0.4975694
0.2403614
0.4711806
0.2018072
0.6607639
0.3090361
第二种做法,用pchip,共52个点,全部有效
首先保存一个M文件:
figure('position',get(0,'screensize'))
axes('position',[0011])
[x,y]=ginput;
然后在commandwindow里输入以下内容:
n=length(x);
s=(1:
n)';
t=(1:
.05:
n)';
u=pchip(s,x,t);
v=pchip(s,y,t);
clfreset
plot(x,y,'.',u,v,'-');
对应的x、y值:
0.5190972
0.8487952
0.5052083
0.7512048
0.4947917
0.6789157
0.5100694
0.6692771
0.5399306
0.7355422
0.5753472
0.8174699
0.596875
0.8620482
0.6190972
0.8777108
0.6149306
0.8138554
0.5878472
0.7427711
0.5878472
0.7427711
0.5635417
0.6716867
0.5350694
0.603012
0.528125
0.563253
0.528125
0.5259036
0.565625
0.5801205
0.6052083
0.6271084
0.634375
0.6186747
0.6190972
0.5716867
0.5878472
0.523494
0.5364583
0.4126506
0.4961806
0.3210843
0.459375
0.2753012
我更喜欢第一种,用spline的,这个能将之间画出弧度,而pchip更像是直接用线段将点依次连接得到的。
使用的是splinetx。
解:
(a)
首先保存一个M文件:
n=3;
xishu=2/(n-1);
x=-1:
xishu:
1;
y=1./(1+25.*x.*x);
fork=1:
n
fenmu=1;
p=1;
forj=1:
n
if(j~=k)
fenmu=fenmu*(x(k)-x(j));
p=conv(p,poly(x(j)));
end
end
c(k,:
)=p*y(k)/fenmu;
end
a=zeros(1,n);
fori=1:
n
forj=1:
n
a(i)=a(i)+c(j,i);
end
end
然后在commandwindow里输入以下内容:
plot(x,y,'*');
holdon;
plot(x,y,'*');
holdon;
x1=-1:
0.01:
1;
y1=0;
fori=1:
n
y1=y1.*x1+a(i)
end
y2=1./(1+25.*x1.*x1);
plot(x1,y1,'b');
holdon;
plot(x1,y2,'g');
即有n=3时,图像:
n=10时,图像:
n=100时,图像:
n=1000时,图像:
可以看出,将[-1,1]做n-1等分的n个插值点,在[-0.92,1]的区间内,随着n趋近于
时
趋近于f(x)。
(b)
先保存M文件:
n=2;
x=2.*rand(n)-1
y=1./(1+25.*x.*x);
n=n^2;
%lagrangc?
?
?
?
?
¡§
fork=1:
n
fenmu=1;
p=1;
forj=1:
n
if(j~=k)
fenmu=fenmu*(x(k)-x(j));
p=conv(p,poly(x(j)));
end
end
c(k,:
)=p*y(k)/fenmu;
end
a=zeros(1,n);
fori=1:
n
forj=1:
n
a(i)=a(i)+c(j,i);
end
end
输出结果:
x=
0.9150-0.6848
0.92980.9412
然后在commandwindow里输入:
plot(x,y,'*');
holdon;
x1=-1:
0.01:
1;
y1=0;
fori=1:
n
y1=y1.*x1+a(i)
end
y2=1./(1+25.*x1.*x1);
plot(x1,y1,'b');
holdon;
plot(x1,y2,'g');
得到以下几幅图:
n=3时,
x=
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趋近于
时
趋近于f(x)。
解:
(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.00040.00000.00000.00000.00000.0000
0.71710.00040.00000.00000.00000.00000.0000
0.75510.00040.00000.00000.00000.00000.0000
0.79500.00040.00000.00000.00000.00000.0000
0.83670.00040.00000.00000.00000.00000.0000
0.88040.00040.00000.00000.00000.00000.0000
0.92610.00050.00000.00000.00000.00000.0000
0.97390.00050.00000.00000.00000.00000.0000
1.02400.00050.00000.00000.00000.00000.0000
Columns8through11
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
0.00000.00000.00000.0000
从数据中可以看出,第3列至第11列均为0,后续的计算无法进行。
(b)
这个复选框是问我们想要哪种的拟合。
(c)
根据多项式来看,P(s)的系数是P(t)的系数的
倍。
没有影响。
t=1900:
10:
2000;
mu=mean(t);
sigma=std(t);
s=(t-mu)./sigma
mu=
1950
sigma=
33.1662
s=
Columns1through6
-1.5076-1.2060-0.9045-0.6030-0.30150
Columns7through11
0.30150.60300.90451.20601.5076
没有更好的取值了。