数值分析作业三次样条插值.docx

数值分析作业三次样条插值.docx

《数值分析作业三次样条插值.docx》由会员分享，可在线阅读，更多相关《数值分析作业三次样条插值.docx（20页珍藏版）》请在冰豆网上搜索。

数值分析作业三次样条插值

数值计算方法作业

实验名称

实验4.3三次样条插值函数

（P126）

4.5三次样条插值函数的收敛性

（P127）

实验时间

姓名

班级

学号

成绩

实验4.3三次样条差值函数

实验目的：

掌握三次样条插值函数的三弯矩方法

实验函数：

2

e2dt

t

X

0.0

0.1

0.2

0.3

0.4

F（X）

0.5000

0.5398

0.5793

「0.6179

0.7554

f（x）=

X

求f（0.13）和f（0.36）的近似值

实验内容:

（1）编程实现求三次样条插值函数的算法，分别考虑不同的边界条件；

（2）计算各插值节点的弯矩值；

（3）在同一坐标系中绘制函数f（x）,插值多项式，三次样条插值多项式的曲线比较插值结果。

实验4.5三次样条差值函数的收敛性

实验目的：

多项式插值不一定是收敛的，即插值的节点多，效果不一定好。

对三次样条插值函数如何呢？

理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。

实验内容：

按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。

实验要求：

（1）随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；

（2）三次样条插值函数的思想最早产生于工业部门。

作为工业应用的例子，考

算法描述:

拉格朗日插值：

错误！

未找到引用源。

nfχ_X）其中错误！

未找到引用源。

是拉格朗日基函数，其表达式为：

Ii（X）L

j=0（Xi一Xj）

牛顿插值：

yixhiMi4hiMiyi

匕一丁）（〒计）,χgxj

式中Mi=S（Xi）.

—hi

hi+hy

-:

hy

'i

hi+hγ

di=-6hGy^比-yhyr=6f[Xij,Xi,Xi.1]hi+hi屮hi屮hi

则Mi满足如下n-1个方程:

.-iMij2M「：

iiMi1=di,i=1,2,...n-1

常用的边界条件有如下几类：

（1）给定区间两端点的斜率mo,m∏,即S（Xo）=y0=mo,S（xn）=yn=mn

（2）给定区间两端点的二阶导数Mo,Mn,即S（XoHy0=M0,S”（Xn）=Mn

（3）假设y=f（X）是以b-a为周期的周期函数，则要求三次样条插值函数S（X）也为周期函数，对S（x）加上周期条件S（P）（XO•0）=S（P）（Xn-0）,p=0,1,2

6Y1—Y0

2M0+M1=—（2j-m0）

对于第一类边界条件有

*h1hi

Mnj2Mn6（mn-也丛）hnhn

2M0+人0M1—d0

FnMnJ+2Mn=dn

Y0=Yn,M°=Mn,S（x°0）=S（Xn-0）,由此推得

'2.7,0......J

M01

d1

田.2..剳.....

M1

d1

M2

d2

....-4..2.Λn4

=

巴..2

MnJ

dn4

<■J

Mn-

I

dn一

那么解就可以为

对于第三类边界条件,

2M00M1%MnT=d0,其中

hl,J0hn,do6（f[Xθ,Xl]-f[Xn-1,Xn]）,那么解就可以为:

hi■hnh1hnh1-hn

程序代码:

1拉格朗日插值函数

Lang.m

funCtiOnf=lang（X,Y,xi）

%X为已知数据的横坐标

%Y为已知数据的纵坐标

%xi插值点处的横坐标

%f求得的拉格朗日插值多项式的值

n=length（X）;

f=0;

fori=1:

n

l=1;

forj=1:

i-1

l=l.*（xi-X（j））∕（X（i）-X（j））;

end;

forj=i+1:

n

l=l.*（xi-X（j））/（X（i）-X（j））;

end;%拉格朗日基函数

f=f+l*Y（i）;

end

fprintf（'%d∖n',f）

return

2牛顿插值函数

newton.m

funCtionf=newton（X,Y,xi）

%X为已知数据的横坐标

%Y为已知数据的纵坐标

%xi插值点处的横坐标

%f求得的拉格朗日插值多项式的值

n=length（X）;

newt=[X',Y'];

%计算差商表

forj=2:

n

fori=n:

-1:

1

ifi>=j

Y（i）=（Y（i）-Y（i-1））∕（X（i）-X（i-j+1））;elseY（i）=0;end

end

newt=[newt,Y'];

end

%计算牛顿插值

f=newt（1,2）;

fori=2:

n

z=1;

fork=1:

i-1

Z=（Xi-X（k））*z;

end

f=f+newt（i-1,i）*z;

end

fprintf（'%d∖n',f）

return

3三次样条插值第一类边界条件

ThreCh.m

funCtionS=ThreCh1（X,Y,dyO,dyn,xi）

%X为已知数据的横坐标

%Y为已知数据的纵坐标

%xi插值点处的横坐标

%S求得的三次样条插值函数的值

%dy0左端点处的一阶导数

%dyn右端点处的一阶导数

n=length（X）-1;

d=zeros（n+1,1）;

h=zeros（1,n-1）;

f1=zeros（1,n-1）;

f2=zeros（1,n-2）;

fori=1:

n%求函数的一阶差商

h（i）=X（i+1）-X（i）;

f1（i）=（Y（i+1）-Y（i））∕h（i）;

end

fori=2:

n%求函数的二阶差商

f2（i）=（f1（i）-f1（i-1））∕（X（i+1）-X（i-1））;d（i）=6*f2（i）;

end

d

（1）=6*（f1

（1）-dyθ）∕h

（1）;

d（n+1）=6*（dyn-f1（n-1））∕h（n-1）;%?

赋初值

A=Zeros（n+1,n+1）;

B=ZeroS（1,n-1）;

C=Zeros（1,n-1）;

fori=1:

n-1

B（i）=h（i）∕（h（i）+h（i+1））;

C（i）=1-B（i）;

end

A（1,2）=1;

A（n+1,n）=1;

fori=1:

n+1

A（i,i）=2;

end

fori=2:

n

A（i,i-1）=B（i-1）;

A（i,i+1）=C（i-1）;

end

M=A∖d;

SymSX;

fori=1:

n

Sx（i）=collect（Y（i）+（f1（i）-（M（i）∕3+M（i+1）∕6）*h（i））*（x-X（i））

+M（i）∕2*（x-X（i））^2+（M（i+1）-M（i））∕（6*h（i））*（x-X（i））^3）;

digits⑷；

Sx（i）=vpa（Sx（i））;%三样条插值函数表达式

end

fori=1:

n

disp（'S（x）='）;

fprintf（'%s（%d,%d）\n',char（Sx（i））,X（i）,X（i+1））;

end

fori=1:

n

ifxi>=X（i）&&xi<=X（i+1）

S=Y（i）+（f1（i）-（M（i）∕3+M（i+1）∕6）*h（i））*（xi-X（i））+M（i）∕2*（xi-X（i））^2+（M（i+1）-M（i））∕（6*h（i））*（xi-X（i））^3;

end

end

disp（'xiS'）;

fprintf（'%d,%d∖n',xi,S）;

return

4三次样条插值第二类边界条件

ThreCh2.m

funCtion[Sx]=Threch2（X,Y,d2yO,d2yn,xi）

X为已知数据的横坐标

%Y为已知数据的纵坐标

%xi插值点处的横坐标

%S求得的三次样条插值函数的值

%d2yθ左端点处的二阶导数

%d2yn右端点处的二阶导数

n=length（X）-1;

d=zeros（n+1,1）;

h=zeros（1,n-1）;

f1=zeros（1,n-1）;

f2=zeros（1,n-2）;

fori=1:

n%求一阶差商

h（i）=X（i+1）-X（i）;

f1（i）=（Y（i+1）-Y（i））∕h（i）;

end

fori=2:

n%求二阶差商

f2（i）=（f1（i）-f1（i-1））∕（X（i+1）-X（i-1））;

d（i）=6*f2（i）;

end

d

（1）=2*d2yθ;

d（n+1）=2*d2yn;%赋初值

A=Zeros（n+1,n+1）;

B=Zeros（1,n-1）;

C=Zeros（1,n-1）;

fori=1:

n-1

B（i）=h（i）∕（h（i）+h（i+1））;

C（i）=1-B（i）;

end

A（1,2）=0;

A（n+1,n）=0;

fori=1:

n+1

A（i,i）=2;

end

fori=2:

n

A（i,i-1）=B（i-1）;

A（i,i+1）=C（i-1）;

end

M=A∖d;

SymSX;

fori=1:

n

Sx（i）=collect（Y（i）+（f1（i）-（M（i）∕3+M（i+1）∕6）*h（i））*（x-X（i））

+M（i）∕2*（x-X（i））^2+（M（i+1）-M（i））∕（6*h（i））*（x-X（i））^3）;

digits⑷；

Sx（i）=vpa（Sx（i））;

end

fori=1:

n

disp（'S（x）='）;

fprintf（'%s（%d,%d）∖n',char（Sx（i））,X（i）,X（i+1））;

end

fori=1:

n

ifxi>=X（i）&&xi<=X（i+1）

S（i）=Y（i）+（f1（i）-（M（i）∕3+M（i+1）∕6）*h（i））*（xi-X（i））+M（i）∕2*（xi-X（i））^2+（M（i+1）-M（i））∕（6*h（i））*（xi-X（i））^3;

end

end

disp（'xiS'）;

fprintf（'%d,%d∖n',xi,S）;

return

5插值节点处的插值结果

main3.m

clear

clc

X=[0.0,0.1,0.2,0.3,0.4];

Y=[0.5000,0.5398,0.5793,0.6179,0.7554];

xi=0.13;

%xi=0.36;

disp（'xi=0.13'）;

%disp（'xi=0.36'）;

disp（'拉格朗日插值结果’）；

Iang（X,Y,xi）;

disp（'牛顿插值结果’）;

newtOn（X,Y,xi）;

disp（'三次样条第一类边界条件插值结果’）;

ThreCh1（X,Y,0.40,0.36,xi）;%0.4,0.36分别为两端点处的一阶导数

disp（'三次样条第二类边界条件插值结果’）;

ThreCh2（X,Y,0,-0.136,xi）;%0,-0.136分别为两端点处的二阶导数

6将多种插值函数即原函数图像画在同一张图上

main2.m

clear

clc

X=[0.0,0.1,0.2,0.3,0.4];

Y=[0.5000,0.5398,0.5793,0.6179,0.7554];

a=linspace（0,0.4,21）;

NUM=21;

L=Zeros（1,NUM）;

N=Zeros（1,NUM）;

S=Zeros（1,NUM）;

B=Zeros（1,NUM）;

fori=1:

NUM

xi=a（i）;

%拉格朗日插值

%牛顿插值

%原函数

%三次样条函数第一类边界条件

P∣ot（a,B,

'--r'

）；

holdon;plot（a,L,

'b'

）；

holdon;

plot（a,N,

holdon;

'r'

）；

'r+'

plot（a,S,

holdon;

legend（'原函数’，’拉格朗日插值’，’牛顿插值’,'三次样条插值’,2）;holdoff

）;

7增加插值节点观察误差变化

main4.m

clear;

clc;

N=5;

%4.5第一问

Ini=zeros（1,1001）;

a=linspace（-1,1,1001）;

Ini=1.∕（1+25*a.^2）;

%节点数量变化次数

fori=1:

3

N=2*N;

t=linspace（-1,1,N+1）;

ft=1.心+25*t^2）;

VaI=IinSPaCe（-1,1,101）;

forj=1:

101

L（j）=lang（t,ft,val（j））;

%插值节点

%插值节点函数值

Sa）=ThreCh1（t,ft,0.074,-0.074,val（j））;end

%三样条第一类边界条件插值

plot（a,Ini,holdon

'k'）%原函数图象

Plot（Val,L,'r'）%拉格朗日插值函数图像

holdon

plot（val,S,'b'）%三次样条插值函数图像

Str=SPrintf（'插值节点为％d时的插值效果',N）;

title（str）;

legend（'原函数’，’拉格朗日插值’，’三次样条插值’）；％显示图例

holdoff

figure

end

8车门曲线

main5.m

clear

clc

X=[0,1,2,3,4,5,6,7,8,9,10];

Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29];

dy0=0.8;

dyn=0.2;

n=length（X）-1;

d=zeros（n+1,1）;

h=zeros（1,n-1）;

f1=zeros（1,n-1）;

f2=zeros（1,n-2）;

fori=1:

nh（i）=X（i+1）-X（i）;

f1（i）=（Y（i+1）-Y（i））∕h（i）;

end

fori=2:

nf2（i）=（f1（i）-f1（i-1））∕（X（i+1）-X（i-1））;

d（i）=6*f2（i）;

end

d

（1）=6*（f1

（1）-dy0）∕h

（1）;

d（n+1）=6*（dyn-f1（n-1））∕h（n-1）;A=Zeros（n+1,n+1）;

B=Zeros（1,n-1）;

C=Zeros（1,n-1）;

fori=1:

n-1

B（i）=h（i）∕（h（i）+h（i+1））;

C（i）=1-B（i）;

end

A（1,2）=1;

A（n+1,n）=1;

fori=1:

n+1

A（i,i）=2;

end

fori=2:

n

A（i,i-1）=B（i-1）;

A（i,i+1）=C（i-1）;

end

M=A∖d;

X=ZeroS（1,n）;

S=ZeroS（1,n）;

fori=1:

n

x（i）=X（i）+0.5;

S（i）=Y（i）+（f1（i）-（M（i）∕3+M（i+1）∕6）*h（i））*（x（i）-X（i））+M（i）∕2*（x（i）-X（i

））^2+（M（i+1）-M（i））∕（6*h（i））*（x（i）-X（i））^3;

end

plot（X,Y,'k'）;holdon;

plot（x,S,'o'）;

title（'三次样条插值效果图’）；

legend（'已知插值节点’，’三次样条插值’）;

holdoff

实验结果：

4.3

1计算插值节点处的函数值

xi=0.13时

Ki=Oe13

5.S36E03e-DOl

半硕⅛5l≡结果

5,665517e-αθl

三样条第——爭⅛⅛⅛界寿√牛価值蜡果

BS=

.500□+.400α*κ-.3S98*x^2+3.698*κ^3

Ξ=

-514Γ-1O,93*κ^3-ι-4,aiε*z-'2-.-⅛O73e-l≠χ<1”OOOOO□e-0O2”ODOoOOe-OO1）

ΞS=

.1O93-∣-39.6≡*κ^3—2G.39⅛2-t-6.O⅛□*κ<2.□C∣O□QOm-OeItA了*OCQDOOe-OO1）

3=

2.476-4Z.93*h,,3+52.4.6⅛"2-17.til*κ¢3.OOUDOOer-OOIJI⅛.ODOooUe-UO1J

χiS

1.30□000e-OOIF5.S31403e-0Ol

三样条第二⅛⅛边界⅛∏f牛插直半吉果

S

.5Q0□÷.380丘7E4⅛'3¢0,1.OOOoOOs-OO1J

ΞS=

5108-9,O7O+κ^3+3.247*sc2-t-.5S73e-l≠κ<1,OOOOOOe-OOU2,OoOOOOe-OOI>

Ξ⅛>=

.I6fi3-bΓ33,93*x^3-22.55*^c-x2÷Ξ-215*x¢2.□OQ□ClDE-（JEJ1,3.□OOΠ口口亡一口门1、

S⅛>=

L.SO7—ΞS.3⅛⅛κ"3-∣-3Ξ.1⅛*≡≡^Ξ-LI*19⅛κ¢3-OOOOODe-OoIJ4-OOOOOOe-OD1>

XlS

1.30000Oe-OOlJllO

5.5□□212—口Olj.>>

Xi=0.36时

拉搭翎曰栖值⅛⅛弟

6,ΓΓ9SS2e-O□I

牛顿⅞G値⅛B耒

T-1633S3θ-O□1

三样毎第一粪边界号lfit⅛⅛値细果

.5□□□-f-*⅛ODO*κ-.3B9B*se^Z-l-3.09S*h□tDj1*OOOOOoe-OOl>

≡M=

“S147—LO.ΘS*3b:

"S-H1*□L3*3T^2—.407Se-L*κ（1*OOOOODe-OO1,2-OOODooe-OOL>

S<«>=

.l□93-∣-39.60≠≥c*3-26.39⅛2+6.O40≠x<2.0OOOOOer-OO1,3.OOOOOOer-OO1>

2.475-⅛T.93⅛3+52.46*x2-17.61*κ¢3T□□OOOOC-OOIJ4,0□0000t-00I）

ZKΞLS

号”0O□QOO⅛-□□Ijβτθβ-i2□O⅛-□□I

三彳羊朵弟二粪边界祭俏插于亘细乐

-r>□OO-I-,38∩5*ZH-1.75^*κ^3

“5103—9.O7r□⅛ιc^3-∣-3x2^2+_5575e-1（1.□OOOOOe-OOIJ2-OOOOoOe-OO1>

“1663-1-33.S3⅛"3—22.SS*3tC2-∣-5.21S*a∑¢2.0□0000⅛-0□lj3.OoOOOo&-OO1）

1.807-2S.84≠x^3-∣-3Ξ.14*x^Ξ-l1.19*x<3.00□000e-0□lp4.00□□□0e-001）

≈ciS

3,60000Oe-OO1,O

U.Ll

e.θ15-45S∈-O□Ij>>I

2将多种插值函数即原函数图像画在同一张图上

0.85

原國数

拉格朗日摘值

牛顿插值摘值

+三样条插值

4.5.1增加插值节点观察误差变化

JSH宜卞熄致为WHT旳JiIn宜孜果

插值节点数为20时的插值效果

从上面三张图可以看出增加插值节点并不能改善差之效果

4.5.2车门曲线