eno matlab,三次样条插值matlab程序 含多种边界条件

【实例简介】

自编的三次样条插值matlab程序 含多种边界条件

0.01087

0.14489

0.368

0.17405

0.4485

0.393

0.2878

0.25891

2.1153

1.5294

2.8188

69.141

0.56548

-21.035

73.614

12.794

58.247

512.32

-28.949

-259.9

42279

25

1.5

0.5

8.5

9.5

10

10.5

图表一、三次样条曲线

(2)坐标軸逆时针旋转45度,相当于节点顺时针旋转45度。设(x,y)为旋转前的

坐标,(x,y)为旋转后的坐标,则有

cos0 -sin0x

sing cos

故旋转后的节点坐标为:

theta=-pi/4;

for i=l: length(x1)

x2()=cos(theta)*x1(i)-sin(theta)*y1(i)i

y2(i)=sin(theta)*x1(i)+cos( theta)*y1(-)

end

fprintf('\t\t\tx2\t\tty 2\n')

disp([x2])i

x2

y 2

5.8

5.6905

6.0097

-5.8697

6.562

-6.166

7.1312

6.2826

7.2889

-6.2876

7.8906

-6.1935

8.4612

-5.9157

8.7519

5.6731

端点处的一阶导数为:

v1=(ul+tan(theta))/(1-ul*tan(theta))i

vn=(un+tan(theta))/(1-un*tan(theta))

fprintf('\t\t\tv1\t\t\tvn\n')

disp([vl vn])i

-0.97849

0.9802

则旋转后的三次样条的系数及图像为:

xx2=[x2(1):0.001:x2(end)]

[Yy2 b2 c2 d2]=spline(x2, y2, xx2,1,v1 vn)

fprintf('\t\t\tb2\t\t\tc2\t\t\td2 \n')i

disp([b2c2(1:end-1,1)d2]);

pLot(x2,y2,*b',xx2,yy2,1-,k");

grid on;

b2

c2

d2

-0.97849

0.67221

0.38277

0.74704

0.43138

0.090754

0.35362

0.28102

0.034909

-0.067629

0.22141

0.053338

0.0061747

0.24664

0.0046897

0.3081

0.2551

0.10233

0.6992

0.43028

0.12195

57--+-

-6.2

-64

5.5

6.5

7.5

图表二、旋转后的三次样条曲线

(3)将第(1)题中所得的样条曲线整体旋转-45度(即顺时针旋转45度),并与

第(2)题的曲线画在同一幅图上,得

for i=1: length(xx1)

xx3(i)=cos(theta)*xx1(i)-sin( theta)*yy1(i)

eno yy3()=sin(theta)xxl(i)+cos(theta)*yy1(i)

plot (xx3, yy3,--, xx2, yY2,-,x2, y2,ok ')igrid on;

legend('旋转前','旋转后');

旋转前

旋转后

-1--

-6.4

-6.6

-6.8

95

图表三、旋转前后样条曲线几何比较

比较上图中的两条曲线,易知曲线不重合,故有以下结论

三次样条函数插值不具备几何不变性!

附录

三次样条插值函数 matlab的m文件程序

function lyy bc d]=spline(x, y, xx, flag, vl, vr)

%三次样条插值函数

%(xy)为插值节点,Xx为插值点;

%flag表端点边界条件类型

%且lag=0:自然样条(端点二阶导数为0);

%flag=1:第一类边界条件(端点一阶导数给定);

%红lag=2:第二类边界条件(端点二阶导数给定);

%Vlv表左右端点处的在边界条件值。

样条函数为:Si(x)=yi+bi(xxi)+ci(xxi)^2+di(xxi)^3

%b,Cd分别为各子区间上的系数值

%yy表插值点处的函数值

if length (x)=length(y

error(输入数据应成对!);

end

length(x);

a=zeros(n-1, 1 );

b

a' d=a dx=a

ay

A

(n); B=zeros(n, 1);

for i=1:n-1

i=y(i

dx(i=x(i+1)-xi):

dy(i)=y(i+1)-y(i);

eno

for i=2: n-1

Ai-1)=dx(i-1)

A1i)=2“(dx(i-1)+dx(i);

A(ii+1)=dx(i)

B(i1)=3(dy()/dx()-dy(i-1)/dx(i-1)

%自然样条端点条件(端点二阶导教为零)%

if flag==0

A(1,1)=1

A(n, n)=1

en

0

%端点一阶导数条件%

a

A(1,1)=2*dx(1);A(12)=dx(1)

A(n-1)=dx(n-1A(n,n)=2dx(n-1)

B(11)=3(dy(1dx(1)v1)

B(n, 1)=3*(vr-dy(n-1)dx(n-1);

0/0

%端点二阶导数条件%

if flag==2

A(11)=2;A(n,n)=2;

B(1,1)=vlB(n,1)

en

AB

for i=1:n-1

d(i)=(c(i+1)-c(/(3°dx(i)

b(1)=dy(i)/dx(i)-dx(i)(2“c(i)+c(i+1)

enc

mm nn=size(xx;

yy=zeros(mm, nn);

for i=l: mm'nn

for ii=1: n-1

xx(i)>=x(ii)&xx(i

bre

ak.

yy(i)=a()+b(i)°(xx(i)-×(j)tc()“(xx(i)x()~2+d(j).(xx(i)-x(j)^3;

nd

en

【实例截图】

【核心代码】