三次样条插值作业题.docx

三次样条插值作业题.docx

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

三次样条插值作业题

例1设

为定义在[0,3]上的函数，有下列函数值表：

xi

0

1

2

3

yi

0

0.5

2

1.5

且

，

，试求区间[0,3]上满足上述条件的三次样条插值函数

本算法求解出的三次样条插值函数将写成三弯矩方程的形式：

其中，方程中的系数

，

，

，

将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。

以下为Matlab代码：

%=============================

%本段代码解决作业题的例1

%=============================

clearall

clc

%自变量x与因变量y，两个边界条件的取值

IndVar=[0,1,2,3];

DepVar=[0,0.5,2,1.5];

LeftBoun=0.2;

RightBoun=-1;

%区间长度向量，其各元素为自变量各段的长度

h=zeros（1,length（IndVar）-1）;

fori=1:

length（IndVar）-1

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

end

%为向量μ赋值

mu=zeros（1,length（h））;

fori=1:

length（mu）-1

mu（i）=h（i）/（h（i）+h（i+1））;

end

mu（i+1）=1;

%为向量λ赋值

lambda=zeros（1,length（h））;

lambda

（1）=1;

fori=2:

length（lambda）

lambda（i）=h（i）/（h（i-1）+h（i））;

end

%为向量d赋值

d=zeros（1,length（h）+1）;

d

（1）=6*（（DepVar

（2）-DepVar

（1））/（IndVar

（2）-IndVar

（1））-LeftBoun）/h

（1）;

fori=2:

length（h）

a=（DepVar（i）-DepVar（i-1））/（IndVar（i）-IndVar（i-1））;

b=（DepVar（i+1）-DepVar（i））/（IndVar（i+1）-IndVar（i））;

c=（b-a）/（IndVar（i+1）-IndVar（i-1））;

d（i）=6*c;

end

d（i+1）=6*（RightBoun-（DepVar（i+1）-DepVar（i））/（IndVar（i+1）-IndVar（i）））/h（i）;

%为矩阵A赋值

%将主对角线上的元素全部置为2

A=zeros（length（d）,length（d））;

fori=1:

length（d）

A（i,i）=2;

end

%将向量λ的各元素赋给主对角线右侧第一条对角线

fori=1:

length（d）-1

A（i,i+1）=lambda（i）;

end

%将向量d的各元素赋给主对角线左侧第一条对角线

fori=1:

length（d）-1

A（i+1,i）=mu（i）;

end

%求解向量M

M=A\d';

%求解每一段曲线的函数表达式

fori=1:

length（h）

Coefs_1=M（i）/（6*h（i））;

Part_1=conv（Coefs_1,...

conv（[-1,IndVar（i+1）],...

conv（[-1,IndVar（i+1）],[-1,IndVar（i+1）]）））;

S_1=polyval（Part_1,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_2=M（i+1）/（6*h（i））;

Part_2=conv（Coefs_2,...

conv（[1,-IndVar（i）],...

conv（[1,-IndVar（i）],[1,-IndVar（i）]）））;

S_2=polyval（Part_2,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_3=（DepVar（i）-M（i）*h（i）^2/6）/h（i）;

Part_3=conv（Coefs_3,[-1,IndVar（i+1）]）;

S_3=polyval（Part_3,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_4=（DepVar（i+1）-M（i+1）*h（i）^2/6）/h（i）;

Part_4=conv（Coefs_4,[1,-IndVar（i）]）;

S_4=polyval（Part_4,[IndVar（i）:

0.01:

IndVar（i+1）]）;

S=S_1+S_2+S_3+S_4;

plot（[IndVar（i）:

0.01:

IndVar（i+1）],S,'LineWidth',1.25）

%在样条插值曲线的相应位置标注该段曲线的函数表达式

text（i-1,polyval（Part_1,3）,...

['\itS',num2str（i）,'（x）=',num2str（Coefs_1）,'（',num2str（IndVar（i+1））,'-x）^{3}+',...

num2str（Coefs_2）,'（x-',num2str（IndVar（i））,'）^{3}+',num2str（Coefs_3）,...

'（',num2str（IndVar（i+1））,'-x）+',num2str（Coefs_4）,'（x-',num2str（IndVar（i））,'）'],...

'FontName','TimesNewRoman','FontSize',14）

holdon

end

%过x=1和x=2两个横轴点作垂线%

line（[1,1],[2.5,-0.5],'LineStyle','--'）;

line（[2,2],[2.5,-0.5],'LineStyle','--'）;

%为x轴和y轴添加标注

xlabel（'\itx','FontName','TimesNewRoman',...

'FontSize',14,'FontWeight','bold'）;

ylabel（'\its（x）','FontName','TimesNewRoman',...

'Rotation',0,'FontSize',14,'FontWeight','bold'）;

最终，三次样条插值函数s（x）表达式为：

曲线的图像如图所示：

例2已知函数值表：

xi

1

2

4

5

yi

1

3

4

2

试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数

本算法求解出的三次样条插值函数将写成三弯矩方程的形式：

其中，方程中的系数

，

，

，

将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。

以下为Matlab代码：

%=============================

%本段代码解决作业题的例2

%=============================

clearall

clc

%自变量x与因变量y的取值

IndVar=[1,2,4,5];

DepVar=[1,3,4,2];

%区间长度向量，其各元素为自变量各段的长度

h=zeros（1,length（IndVar）-1）;

fori=1:

length（IndVar）-1

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

end

%为向量μ赋值

mu=zeros（1,length（h））;

fori=1:

length（mu）-1

mu（i）=h（i）/（h（i）+h（i+1））;

end

mu（i+1）=0;

%为向量λ赋值

lambda=zeros（1,length（h））;

lambda

（1）=0;

fori=2:

length（lambda）

lambda（i）=h（i）/（h（i-1）+h（i））;

end

%为向量d赋值

d=zeros（1,length（h）+1）;

d

（1）=0;

fori=2:

length（h）

a=（DepVar（i）-DepVar（i-1））/（IndVar（i）-IndVar（i-1））;

b=（DepVar（i+1）-DepVar（i））/（IndVar（i+1）-IndVar（i））;

c=（b-a）/（IndVar（i+1）-IndVar（i-1））;

d（i）=6*c;

end

d（i+1）=0;

%为矩阵A赋值

%将主对角线上的元素全部置为2

A=zeros（length（d）,length（d））;

fori=1:

length（d）

A（i,i）=2;

end

%将向量λ的各元素赋给主对角线右侧第一条对角线

fori=1:

length（d）-1

A（i,i+1）=lambda（i）;

end

%将向量d的各元素赋给主对角线左侧第一条对角线

fori=1:

length（d）-1

A（i+1,i）=mu（i）;

end

%求解向量M

M=A\d';

%求解每一段曲线的函数表达式

fori=1:

length（h）

Coefs_1=M（i）/（6*h（i））;

Part_1=conv（Coefs_1,...

conv（[-1,IndVar（i+1）],...

conv（[-1,IndVar（i+1）],[-1,IndVar（i+1）]）））;

S_1=polyval（Part_1,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_2=M（i+1）/（6*h（i））;

Part_2=conv（Coefs_2,...

conv（[1,-IndVar（i）],...

conv（[1,-IndVar（i）],[1,-IndVar（i）]）））;

S_2=polyval（Part_2,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_3=（DepVar（i）-M（i）*h（i）^2/6）/h（i）;

Part_3=conv（Coefs_3,[-1,IndVar（i+1）]）;

S_3=polyval（Part_3,[IndVar（i）:

0.01:

IndVar（i+1）]）;

Coefs_4=（DepVar（i+1）-M（i+1）*h（i）^2/6）/h（i）;

Part_4=conv（Coefs_4,[1,-IndVar（i）]）;

S_4=polyval（Part_4,[IndVar（i）:

0.01:

IndVar（i+1）]）;

S=S_1+S_2+S_3+S_4;

plot（[IndVar（i）:

0.01:

IndVar（i+1）],S,'LineWidth',1.25）

%在样条插值曲线的相应位置标注该段曲线的函数表达式

text（i,polyval（Part_1,5）,...

['\itS',num2str（i）,'（x）=',num2str（Coefs_1）,'（',num2str（IndVar（i+1））,'-x）^{3}+',...

num2str（Coefs_2）,'（x-',num2str（IndVar（i））,'）^{3}+',num2str（Coefs_3）,...

'（',num2str（IndVar（i+1））,'-x）+',num2str（Coefs_4）,'（x-',num2str（IndVar（i））,'）'],...

'FontName','TimesNewRoman','FontSize',14）

holdon

end

%过x=2和x=4两个横轴点作垂线%

line（[2,2],[4.5,0.5],'LineStyle','--'）;

line（[4,4],[4.5,0.5],'LineStyle','--'）;

%为x轴和y轴添加标注

xlabel（'\itx','FontName','TimesNewRoman',...

'FontSize',14,'FontWeight','bold'）;

ylabel（'\its（x）','FontName','TimesNewRoman',...

'Rotation',0,'FontSize',14,'FontWeight','bold'）;

最终，三次自然样条插值函数s（x）表达式为：

曲线的图像如图所示：

例3课后习题与思考题第7题

xi

0.25

0.30

0.39

0.45

0.53

yi

0.5000

0.5477

0.6245

0.6708

0.7280

试求在区间[0.25,0.53]上满足上述函数表所给出的插值条件的三次自然样条插值函数s（x）

求解出的三次样条插值函数将写成三弯矩方程的形式：

本题采用和例2基本相同的Matlab代码，只改变初始条件。

最终，三次自然样条插值函数s（x）表达式为：

曲线的图像如图所示：

例4课后习题与思考题第6题

求

的二次插值式

，使：

并计算

的近似值并估计误差。

本题采用拉格朗日二次插值法进行计算：

以下为Matlab代码：

%=============================

%本段代码解决课本第2章习题与思考题第6题

%=============================

clearall

clc

%自变量x与因变量y的取值

IndVar=[100,121,144];

DepVar=[10,11,12];

%构造拉格朗日插值函数

Coefs_1=DepVar

（1）/...

（（IndVar

（1）-IndVar

（2））*（IndVar

（1）-IndVar（3）））;

Part_1=conv（Coefs_1,...

conv（[1,-IndVar

（2）],[1,-IndVar（3）]））;

f_1=polyval（Part_1,[IndVar

（1）:

0.01:

IndVar（3）]）;

Coefs_2=DepVar

（2）/...

（（IndVar

（2）-IndVar

（1））*（IndVar

（2）-IndVar（3）））;

Part_2=conv（Coefs_2,...

conv（[1,-IndVar

（1）],[1,-IndVar（3）]））;

f_2=polyval（Part_2,[IndVar

（1）:

0.01:

IndVar（3）]）;

Coefs_3=DepVar（3）/...

（（IndVar（3）-IndVar

（1））*（IndVar（3）-IndVar

（2）））;

Part_3=conv（Coefs_3,...

conv（[1,-IndVar

（1）],[1,-IndVar

（2）]））;

f_3=polyval（Part_3,[IndVar

（1）:

0.01:

IndVar（3）]）;

f=f_1+f_2+f_3;

plot（[IndVar

（1）:

0.01:

IndVar（3）],f,'LineWidth',1.25）;

%在样条插值曲线的相应位置标注该段曲线的函数表达式

text（110,polyval（Part_1,110）+polyval（Part_2,110）+polyval（Part_3,110）,...

['\itP_2（x）=',num2str（Coefs_1）,'（x-',num2str（IndVar

（2））,...

'）（x-',num2str（IndVar（3））,'）+',num2str（Coefs_2）,'（x-',num2str（IndVar

（1））,...

'）（x-',num2str（IndVar（3））,'）+',num2str（Coefs_3）,'（x-',num2str（IndVar

（1））,...

'）（x-',num2str（IndVar

（2））,'）'],...

'FontName','TimesNewRoman','FontSize',14）

%为x轴和y轴添加标注

xlabel（'\itx','FontName','TimesNewRoman',...

'FontSize',14,'FontWeight','bold'）;

ylabel（'\itP_2（x）','FontName','TimesNewRoman',...

'Rotation',0,'FontSize',14,'FontWeight','bold'）;

%求115平方根的估计值

fun115=polyval（Part_1,115）+polyval（Part_2,115）+polyval（Part_3,115）

%求115平方根的准确值

sqrt（115）

%求误差，单位为%

est=100*（fun115-sqrt（115））/fun115

最终，拉格朗日插值函数

的表达式为：

从计算结果可以看出，

的准确值为10.7238，而通过拉格朗日插值法求出的近似值为10.7228，误差为0.0098%

的图像为：