昆明理工数值分析Ronge现象

课题六 Ronge 现象的产生与克服

一、目的和意义

1、深刻认识多项式插值的优缺点；

2、明确插值方法的不稳定性如何克服；

3、体会精度与节点，插值方法的关系；

4、利用计算机绘图来显示插值函数，使结果形象化；

5、理解三次样条插值在实际应用中的价值。

二、问题提出

给定函数

()2

11x x f +=

，55≤≤-x ，及节点N j N j x j ,,2,1,0,10

5 =⨯+-=，

试用何种插值方法可克服Ronge 现象。 1、用多项式插值计算出下列插值。

()

0y p ()

k p y ⨯+5.045.0

()

k p y ⨯--=5.045.0 9,,2,1 =k ，观察是否会产生

Ronge 现象。

2、用下列插值方法进行计算，并比较它们克服Ronge 现象的效果。 （1）分段线性插值；

（2）三次样条函数插值（一），条件共为四个，任意选择两个。此处选择前两个为：

()()j j y x f x s = N j ,,2,1,0 = ()()00x f x s y '=' ()()y y y x f x s '='

三、计算方法及公式

1、对于多项式插值

采用拉格朗日插值法，n 次插值基函数为 )

())(()()

())(()()(11111111++-++---------=

n k k k k k k n k k k x x x x x x x x x x x x x x x x x l

其中，k=1，1，…，n+1 拉格朗日差值多项式可表示为

∑+==1

1

)()(n k k k n x l y x L

注：由于matlab 数组的首地址都以1开始，故该插值公式也从1开始。 计算步骤：（1）输入插值节点及节点值。 （2）循环计算各次插值基函数值。

（3）累加基函数及对应节点乘积，得到插值结果。 2、分段线性插值

根据分段线性插值的定义，在每个区间[x i ,x i+1]上，插值函数可表示为 11

11()

i i i i i i

x x x x h i i x x x x I x f f +++--+--=

+

计算步骤：（1）根据输入节点及节点函数值，建立各个区间的插值函数。 （2）根据要插值的节点，判断其所在区间，利用该区间的插值函

数计算其插值结果。

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

三次样条插值函数在每个区间[x i ,x i+1]上，插值函数的表达式如下：

33

22

111()()11666

6

()()

()

i i i i i i x x x x M h x x M h x x i i i i h

h

h

h

S x M

M

f f +++----++=++-

+-

其中M i ,M i+1为节点i 及i+1处的弯矩值，为未知量。

M 的计算方法如下：

由三次样条函数导数连续的性质，可有以下关系式：

112i i i i i i M M M d μλ-+++=

i=2,3,…n-1。

由于所提供节点为等距节点，故 1/2i μ= 1/2i λ= 11

2

23

i i i f f f i h d +--+=

由于边界条件为第一类边界条件，故有

21

61212(')f f h h M M f -+=-

1

6

12(')n n f f n n n h h

M M f ---+=-

因此，可以组成三对角方程组利用追赶法对M 进行求解。 计算步骤：（1）根据已知节点值及边界点一阶导数值，组成M 求解的系数矩

阵以及求解向量d 。

（2）利用追赶法求解M 。

（3）根据要插值的节点，判断其所在区间，利用该区间的插值函

数计算其插值结果。

四、程序设计

1、拉格朗日插值程序

function mainfunction()

x=-5:5;

y=1./(1+x.^2);

xx1=-4.95:0.5:-0.45;xx2=0;xx3=0.45:0.5:4.95;

xx=[xx1 xx2 xx3];

yy=lglr(x,y,xx); %调用拉格朗日插值子程序

xx'

yy'

xxx=-5:0.1:5;

yf=1./(1+xxx.^2);

plot(xxx,yf,'k.-',xx,yy)

xlabel('x轴')

ylabel('y轴')

gtext('理论曲线')

gtext('拉格朗日插值曲线')

function yy=lglr(x,y,xx)%拉格朗日插值子程序[n,m]=size(x);[n,k]=size(xx);

yy=zeros(1,k);

for ii=1:k

for i=1:m

L=y(i);

for j=1:m

if i~=j

L=L*(xx(ii)-x(j))/(x(i)-x(j));

end

end

yy(ii)=yy(ii)+L;

end

end

2、分段线性插值程序

function mainfunction()

x=-5:5;

y=1./(1+x.^2);

xx1=-4.95:0.5:-0.45;xx2=0;xx3=0.45:0.5:4.95;

xx=[xx1 xx2 xx3];

yy=fdx(x,y,xx); %调用分段线性插值子程序xxx=-5:0.1:5;

yf=1./(1+xxx.^2);

plot(xxx,yf,'k.-',xx,yy,'kv-')

xlabel('x轴')

ylabel('y轴')