数学实验“等距节点插值,Hermite插值,分段插值(线性,二次,三次)”实验报告(内含matlab程序)

西京学院数学软件实验任务书

实验十六实验报告

一、实验名称：等距节点插值，Hermite 插值，分段插值（线性，二次，三次）。

二、实验目的：进一步熟悉等距节点插值，Hermite 插值，分段插值（线性，二次，三次）。

三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。 四、实验原理：

1．等距节点插值：

差分分为前向差分、后向差分和中心差分三种，它们的记法及定义如下所示:

n 阶前向差分公式111()()()n n n i i i f x f x f x --+?=?-? n 阶后向差分公式111()()()n n n i i i f x f x f x ---?=?-?

n 阶中心差分公式11112

2

()()()n n n i i i f x f x f x δδδ--+-=-

其中：? -前向差分；? -后向差分；δ -中心差分。 假设000()()()()i i i i f x f x f x f x δ?=?==，为了方便计算，构造差分表(()i i f f x =)。

这里只说明前向牛顿插值，其多项式可表示为如下形式：

0()()N x N x th =+

20000()()()()12n t t t f x f x f x f x n

??

??

??

=+?+?++? ? ? ???

??

??

其中h 为步长，10h x x =-，且的取值范围为0t n ≤≤。 2．埃尔米特插值：

埃尔米特插值法满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值，对于有高阶导数的情况，埃尔米特插值多项式比较复杂，在实际应用中，常常遇到的是函数值与一阶导数值给定的情况，在这种情况下，n 个节点12,,n x x x 的埃尔米特插值多项式()H x 的表达形式如下所示：

1()[()(2)]n

i i i i i i i H x h x x a y y y ==--+∑

其中(),''()i i i i y y x y y x == 2

111

(

),n

n

j i i j j i j

i j

j i

j i x x h a x x x x ==≠≠-==--∑

∏

3．分段插值：

给定插值节点i x 、节点函数值i y 及对应的导数值

'(0,1,2,,)i y i N = ，则满足下面条件

(),'()'i i i i p x y p x y ==

的分段埃尔米特插值函数()p x 的表达式如下所示：

22

11122111()(12)()(12)()'()(

)'()()i i i i i i i i i i i

i i i i i i i i

x x x x x x x x p x y y h h h h x x x x y x x y x x h h ++++++----=++++----+-

11,(0,1,2,1),[,]i i i i i h x x i N x x x ++=-=-∈

五、实验内容：

%等距节点插值

function [f,f0]= dengjujiedian(x,y,x0)

syms t ;

if (length(x) == length(y)) n = length(x); c(1:n) = 0.0; else

disp('x 和y 的维数不相等！'); return ; end

f = y(1); y1 = 0;

xx =linspace(x(1),x(n),(x(2)-x(1))); if (xx ~= x)

disp('节点之间不是等距的！'); return ; end

for (i=1:n-1)

for (j=1:n-i)

y1(j) = y(j+1)-y(j); end

c(i) = y1(1); l = t;

for (k=1:i-1)

l = l*(t-k); end ;

f = f + c(i)*l/factorial(i); simplify(f);

y = y1;

end

f0=subs(f,'t',(x0-x(1))/(x(2)-x(1)));

%埃尔米特插值

function [f,f0]= Hermite(x,y,y_1,x0)

syms t;

f = 0.0;

if(length(x) == length(y))

if(length(y) == length(y_1))

n = length(x);

else

disp('y和y的导数的维数不相等！');

return;

end

else

disp('x和y的维数不相等！');

return;

end

for i=1:n

h = 1.0;

a = 0.0;

for j=1:n

if( j ~= i)

h = h*(t-x(j))^2/((x(i)-x(j))^2);

a = a + 1/(x(i)-x(j));

end

end

f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i)); end

f0=subs(f,'t',x0);

%分段差值

function [f,f0] = fenduan(x,y,y_1,x0)

syms t;

f = 0.0;

f0 = 0.0;

if(length(x) == length(y))

if(length(y) == length(y_1))

n = length(x);

else

disp('y和y的导数的维数不相等！');

return;

end

else

disp('x和y的维数不相等！');

return;

end

for i=1:n

if(x(i)<=x0)&& (x(i+1)>=x0)

index = i;

break;

end

end

h = x(index+1) - x(index);

fl = y(index)*(1+2*(t-x(index))/h)*(t-x(index+1))^2/h/h + ...

y(index+1)*(1-2*(t-x(index+1))/h)*(t-x(index))^2/h/h;

fr = y_1(index)*(t-x(index))*(t-x(index+1))^2/h/h + ...

y_1(index+1)*(t-x(index+1))*(t-x(index))^2/h/h;

f = fl + fr;

f0 = subs(f,'t',x0);