计算方法大作业1 克服Runge现象

S ( x)

( xi x)3 6 hi 1
M i1

( x xi1 )3 6 hi 1
Mi


yi 1

h2 i 1 6
M i1

xi x hi 1

Байду номын сангаас
yi

h2 i 1 6
Mi

x
xi1 hi 1
, x [ xi1, xi ]
x3
x2
x
1
S1 ( x)
-0.34685
0.2086
0.073964
0.038462
S2 (x)
因此考虑三种边界条件的三弯矩方程分别如下：
（1）第一边界条件
2 1


2
2
2
M1 d1 1M 0

M2


d2





n2
2
n2


M
n2


dn2


n 1
2


M
n
1

在给定 n+1 个节点和相应的函数值后构造 n 次的 Lagrange 插值多项式，其等距节 点的高次多项式插值函数在端点会产生剧烈震荡，因此多项式对函数的逼近并不是随着 次数增高而越理想，这种现象就是 Runge 现象。解决 Runge 现象的方法通常有分段线性 插值、三次样条插值等方法。这里采用三次样条插值方法。
由式（1）内点拼接条件，可得
i M i1 2M i i M i1 d j i 1, 2,..., n 1
（3） （4）
其中
i

hi 1 hi1
, hi

i

hi hi 1

,d hi
i
6
yi 1 yi hi
yi


h1i
y1i
/h(1i h i ) .
Lagrange 插值多项式公式，对插值区间内任意 x 的函数值 y 可通过下式求得
y(x)
n
n
yk (
k 1
j 1
x xj ) xk x j
jk
采用 MATLAB 编写 Lagrange 函数：
Lagrange.m
function y= lagrange(x0,y0,x) n=length(x0); m=length(x); for i=1:m

d n1

n 1 M
n

（2）第二边界条件
（5）
2
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
（3）第三边界条件
2 1

1
2
1
M0 d0

M1


d1




n 1
2
n
1


M
n
1
z=x(i); s=0.0; for k=1:n
p=1.0; for j=1:n
if j~=k p=p*(z-x0(j))/(x0(k)-x0(j));
end end s=p*y0(k)+s; end y(i)=s; end
为了观察Runge现象，在MATLAB命令窗口中输入（n=10）：
>> x=linspace(-1,1,10); >> y=[1./(1+25*(x.^2))]; >> x0=[-1:0.01:1]; >> y0=lagrange(x,y,x0); >> y1=1./(1+25*(x0.^2)); >> plot(x0,y0,'--r') >> hold on >> plot(x0,y1,'-b')

dn1

1
2


Mn


dn

（6）
2 1


2
2
2
1 M1 d1

M2


d2




n 1
2
n
1


M
n
1

dn1
n
n 2 M n dn
三次样条插值克服 Runge 现象
在命令窗口输入：
clear；clc >> x=linspace(-1,1,11); >> y=[1./(1+25*(x.^2))]; >> dx0=0.073964497; >> dxn=-0.073964497; >> S=scfit(x,y,dx0,dxn); >> x1=-1:0.001:-0.8; y1=polyval(S(1,:),x1-X(1)); >> x2=-0.8:0.001:-0.6; y2=polyval(S(2,:),x2-x(2)); >> x3=-0.6:0.001:-0.4; y3=polyval(S(3,:),x3-x(3)); >> x4=-0.4:0.001:-0.2; y4=polyval(S(4,:),x4-x(4)); >> x5=-0.2:0.001:0; y5=polyval(S(5,:),x5-x(5)); >> x6=0:0.001:0.2; y6=polyval(S(6,:),x6-x(6)); >> x7=0.2:0.001:0.4; y7=polyval(S(7,:),x7-x(7)); >> x8=0.4:0.001:0.6; y8=polyval(S(8,:),x8-x(8)); >> x9=0.6:0.001:0.8; y9=polyval(S(9,:),x9-x(9)); >> x10=0.8:0.001:1.0; y10=polyval(S(10,:),x10-x(10));
五、 计算结果
从图 1 与图 2 可以看出，采用多项式插值法时，其等距节点的高次多项式插值函数
在端点会产生剧烈震荡，且震荡随着 n 的增加会越明显，即 Runge 现象。
图 1 n=10 时 Lagrange 插值曲线
图 2 n=5，10，15 时 Lagrange 插值曲线
为了克服 Runge 现象我们采用了三次样条插值（第二类边界条件），从图 3 可以看
（2）
要得到 S (x) 的表达式，首先要求出 M i 的值，对（2）式求一阶导，得
S '(x)

( xi x)2 2 hi 1
M i1

( x xi1)2 2 hi 1
Mi

yi yi1 hi 1

1 6 (Mi
M i1 )hi1, x [ xi1, xi ]
6
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
出，此时不会出现多项式插值时出现的 Runge 现象，插值函数与原函数拟合程度很高， 插值效果明显提高。
图 3 三次样条插值曲线
经过编程求解可以得到三次样条插值函数 Sn (x) ，最终得到的三次样条插值多项式的 系数如下表。
S ( x)
设在区间[a, b]上给定 n+1 个节点 xi （a = x0 < x1 < … < xn = b），在节点 xi 处的函数
值为 yi f (xi ) (i = 0,1,…,n).若函数 S (x) 满足以下三个条件： （1）在每个子区间[ xi , xi1 ] (i = 0,1,…,n-1)上， S (x) 是三次多项式；
（2） S (xi ) yi (i = 0,1,…,n)； （3）在区间[a, b]上， S (x) 的二阶导数 S ''(x) 连续；
则称 S (x) 为函数 y f (x) 在区间[a, b]上的三次样条插值函数。
由定义知要求出 S (x) ，在每个小区间[ xi , xi1 ] 上要确定 4 个待定系数，而共有 n 个小区间，故要确定 4n 个参数。根据 S (x) 在[a, b]上二阶导数连续，因此有
逼近效果往往不理想，会产生 Runge 现象。采用三次样条插值可以有效克服这一现象， 并且取得较好的插值效果。
采用函数 f (x) 1 进行数值插值，插值区间为[-1,1]，取不同的节点数 n，在 1 25x 2
区间[-1,1]上取等距节点为插值点，给定节点为 xj=-1+jh,h=0.1，j=0,…，n。再采用三次 样条插值，并比较二者插值效果。 二、 数学原理
界条件，常见的有以下三种：
1
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
（1） 已知两端的二阶导数值 S ''(x0 ) f ''0, S ''(xn ) f ''n ； （2） 已知两端的一阶导数值 S '(x0 ) f '0, S '(xn ) f 'n ；
（3） 当 f (x) 是以 xn x0 为周期的周期函数时，则要求 S (x) 也是周期函数。边界
S (xi 0 ) S x(i 0 )

S
'
(xi

0) S
xi' (
0 )i

S
'
'
x(i

0)S
xi' ' (
0)
1 ,n2, . . . , 1
（1）
这里共有了 3n-3 个条件，再加上条件（2）中的 n+1 个插值条件，共有 4n-2 个条件，
因此还需要 2 个方程才能确定 S (x) .通常可在区间[a, b]的端点 a x0,b xn 上各加一个边
temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end
采用三次样条插值，构造函数：
scfit.m
function S=scfit(X,Y,dx0,dxn) N=length(X)-1; H=diff(X);
4
《计算方法》大作业（一）
D=diff(Y)./H; A=H(2:N-1); B=2*(H(1,N-1)+H(2:N)); C=H(2:N); U=6*diff(D); B(1)=B(1)-H(1)/2; U(1)=U(1)-3*(D(1)); B(N-1)=B(N-1)-H(N)/2; U(N-1)=U(N-1)-3*(-D(N)); for k=2:N-1
条件满足
S (x0 0) S (xn

S
''( x0

0)

S
''(
xn
0), S 0)
'( x0

0)

S
'( xn

0)
.
设 S (x) 在 xi 处的二阶导数值为 M i (i = 0,1,…,n)， M i 为待定参数，由式（1）可以得
到 S (x) 在区间[ xi1 , xi ] 上的表达式为
（7）
三次样条插值函数的关键是要解式（5）、（6）或（7）所示的方程组，将得到的M i
代入公式（2）中可得到相应的插值函数。 三、 算法框图
输入插值节点个数 n 插值 节点 ,插值点函数值
第一类 输入
计算 选择边界条件
第几类边界条件 第三类
第二类 输入
得到 n-1 阶方程组 式（5）
得到 n 阶方程组 式（7）
5
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
>> x=-1:0.001:1; >> y=1./(25*x.^2+1); >> plot(x,y,'r') >> hold on >>plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,x10,y10,x,y,'-')
用追赶法 求解式（ 5）（6）LU 分解法求解式（7）得到 Mi
将 Mi 代入式（2）得到三次样 条插值函数 3
得到 n+1 阶方程组 式（6）
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
四、 程序实现
对给定的 n 个插值节点 x1, x2 ,..., xn 及对应的函数值 y1, y2 ,..., yn ，利用（n-1）次
《计算方法》
上机大作业（一）
题目：三次样条插值法克服 Runge 现象
姓名： 李楠 班级： B1603193 学号： 116031910075
2016 年 11 月
《计算方法》大作业（一）
三次样条插值克服 Runge 现象
一、 题目描述 对于多项式插值，插值多项式的次数随着节点个数的增加而升高，然而高次插值的