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

-3.0000 15.0000 -18.0000 0
0.5000 -2.0000 1.5000 0
2.6667 -8.0000 5.3333 0
a =
0.8333 0 0 0
a =
-2.1667 0 0 0
a =
-1.6667 0 0 0
a =
1 0 0 0
a =
1 -5 0 0
a =
1 10 0 0
x1=-1:0.01:1;
y1=0;
for i=1:n
y1=y1.*x1+a(i)
end
y2=1./(1+25.*x1.*x1);
plot(x1,y1,'b');
hold on;
plot(x1,y2,'g');
得到以下几幅图：
n=3时，
x =
0.4121 -0.9077 0.3897
-0.9363 -0.8057 -0.3658
0.8407 0.4022 0.1224 -0.2475 -0.4191 0.7593 -0.8115
-0.8946 0.3327 0.7637 -0.6182 0.2342 0.6355 0.1970
0.4757 0.0783 0.3384 -0.1435 -0.4694 -0.4785 -0.0582
axes('position',[0 0 1 1])
[x,y] = ginput;
然后在command window里输入以下内容：
n = length(x);
s = (1:n)';
t = (1:.05:n)';
u = spline(s,x,t);
v = spline(s,y,t);
clf reset
-0.4462 0.6469 0.9004
n=10时，
x =
Columns 1 through 7
-0.3210 0.9661 -0.6578 0.7110 0.1660 -0.7845 -0.6425
0.9033 -0.3971 -0.9348 0.2895 -0.4964 0.8126 -0.1542
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更像是直接用线段将点依次连接得到的。
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
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'))
(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)));
Baidu Nhomakorabeaend
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
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.9298 0.9412
然后在command window里输入：
plot(x,y,'*');
hold on;
a =
1 8 0 0
a =
1 0 0 0
a =
1.0000 0 9.1667 0
a =
1.0000 0 -8.8333 0
a =
1.0000 0 -7.3333 0
a =
1.0000 0 -2.0000 0
a =
1.0000 0 -2.0000 -5.0000
a =
1.0000 0 -2.0000 -5.0000
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
使用的是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)));
1
p =
1 0
fenmu =
-1
p =
1 -2 0
fenmu =
2
p =
1 -5 6 0
c =
0.8333 -5.0000 9.1667 -5.0000
-3.0000 15.0000 -18.0000 0
fenmu =
2
p =
1 0
fenmu =
2
p =
1 -1 0
fenmu =
-2
p =
1 -4 3 0
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.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.7887 -0.0852 0.3903
n=10时，数据太大，没运行出来。
可以看出，将[-1,1]做n-1等分的n个插值点，在[-0.92,0.92]的区间内，随着n趋近于 时 趋近于f(x)。
解：
(a)
t = 1900:10:2000
V = vander(t)
输出结果：
t =
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 -3 2
fenmu =
-6
p =
1 -6 11 -6
c =
0.8333 -5.0000 9.1667 -5.0000
fenmu =
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
x1=-1:0.01:1;
y1=0;
for i=1:n
y1=y1.*x1+a(i)
end
y2=1./(1+25.*x1.*x1);
plot(x1,y1,'b');
hold on;
plot(x1,y2,'g');
即有n=3时，图像：
n=10时，图像：
n=100时，图像：
n=1000时，图像：
可以看出，将[-1,1]做n-1等分的n个插值点，在[-0.92,1]的区间内，随着n趋近于 时 趋近于f(x)。
c =
0.8333 -5.0000 9.1667 -5.0000
-3.0000 15.0000 -18.0000 0
0.5000 -2.0000 1.5000 0
fenmu =
3
p =
1 0
fenmu =
6
p =
1 -1 0
fenmu =
6
p =
1 -3 2 0
c =
0.8333 -5.0000 9.1667 -5.0000
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.8855 -0.7440 0.9633 -0.5476 -0.3122 -0.3746 -0.9328
-0.1645 0.9982 -0.6872 -0.2308 0.1681 -0.6770 -0.8624
Columns 8 through 10
-0.3608 0.2219 0.7507
0.0617 0.5576 0.0361
axes('position',[0 0 1 1])
[x,y] = ginput;
然后在command window里输入以下内容：
n = length(x);
s = (1:n)';
t = (1:.05:n)';
u = pchip (s,x,t);
v = pchip (s,y,t);
clf reset
-0.4618 0.3962 -0.6191 -0.0360 0.6488 0.1887 0.3919
-0.1543 0.3331 -0.2622 -0.7588 0.9653 -0.9550 0.3998
0.0957 -0.6437 -0.0785 0.1790 0.4605 -0.1495 0.2771
0.3089 -0.1531 0.8872
-0.1848 -0.8184 0.2754
0.6400 -0.4671 0.9154
0.4367 -0.6927 -0.5186
0.9373 -0.4380 0.3522
0.0627 -0.1198 -0.4219
-0.3497 0.0543 0.3436
k=j;
end
end
在command window中输入：
s=u-x(k);
v=y(k)+s.*delta(k)
输出结果：
v =
15.6000
解：
第一种做法，用spline，共55个点，其中，54个有效
首先保存你一个M文件：
figure('position',get(0,'screensize'))
a =
1.0000 0 -2.0000 -5.0000
a =
1.0000 0 -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)<u
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
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
然后在command window里输入以下内容：
plot(x,y,'*');
hold on;
plot(x,y,'*');
hold on;
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