问题背景
在数值分析这个数学分支中,样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差,这样就避免了使用高阶多项式所出现的龙格现象,所以样条插值得到了流行。 该问题的方法很好地解决了第一问的龙格现象,得到的插值函数有着很好的性质。 样条的概念来自船体放样中中所用的木头长条。一条船的外曲线要求很光顺,又要通过一些点, 这些点称为型值点,工人技术人员就用这种样条来画曲线。两个相邻型值点之间的曲线,近似地为一个三次多项式。整个曲线通过各型值点处,其一阶导数、二阶导数是连续的。
数值演算
由matlab编写程序,利用D1样条插值函数在4到6之间的多项式表达式求值,得到的结果为:1.609770287689208 用D1样条插值得到的曲线为:
源代码
clear;clc
;
format long
x
=[1 2 3 4 6
];
f
=[log
(1
),log
(2
),log
(3
),log
(4
),log
(6
)];
f11
=1
;
f55
=1/6
;
f12
=(f
(2
)-f
(1
))/
(x
(2
)-x
(1
));
f23
=(f
(3
)-f
(2
))/
(x
(3
)-x
(2
));
f34
=(f
(4
)-f
(3
))/
(x
(4
)-x
(3
));
f45
=(f
(5
)-f
(4
))/
(x
(5
)-x
(4
));
f123
=(f23-f12
)/
(x
(3
)-x
(1
));
f234
=(f34-f23
)/
(x
(4
)-x
(2
));
f345
=(f45-f34
)/
(x
(5
)-x
(3
));
f112
=(f11-f12
)/
(x
(1
)-x
(2
));
f455
=(f55-f45
)/
(x
(5
)-x
(4
));
for i
=2:4
u
(i
)=(x
(i
)-x
(i-1
))/
(x
(i+1
)-x
(i-1
));
t
(i
)=(x
(i+1
)-x
(i
))/
(x
(i+1
)-x
(i-1
));
end
A
=[2 1 0 0 0
;u
(2
) 2 t
(2
) 0 0
;0 u
(3
) 2 t
(3
) 0
;0 0 u
(4
) 2 t
(4
);0 0 0 1 2
];
b
=6*
[f112
;f123
;f234
;f345
;f455
];
M
=A\b
;
s5
=f
(4
)+
(f
(5
)-f
(4
))/
(x
(5
)-x
(4
))*
(5-x
(4
))-
(M
(5
)/6+M
(4
)/3
)*
(x
(5
)-x
(4
))*
(5-x
(4
))+M
(4
)/2*
(5-x
(4
))^2+
(M
(5
)-M
(4
))/
(x
(5
)-x
(4
))*
(5-x
(4
))^3/6
X
=linspace
(1,2,100
);
Y
=f
(1
)+
(f
(2
)-f
(1
))/
(x
(2
)-x
(1
)).*
(X-x
(1
))-
(M
(2
)/6+M
(1
)/3
)*
(x
(2
)-x
(1
)).*
(X-x
(1
))+M
(1
)/2.*
(X-x
(1
)).^2+
(M
(2
)-M
(1
))/
(x
(2
)-x
(1
)).*
(X-x
(1
)).^3/6
;
plot
(X,Y
);
hold on
;
X
=linspace
(2,3,100
);
Y
=f
(2
)+
(f
(3
)-f
(2
))/
(x
(3
)-x
(2
)).*
(X-x
(2
))-
(M
(3
)/6+M
(2
)/3
)*
(x
(3
)-x
(2
)).*
(X-x
(2
))+M
(2
)/2.*
(X-x
(2
)).^2+
(M
(3
)-M
(2
))/
(x
(3
)-x
(2
)).*
(X-x
(2
)).^3/6
;
plot
(X,Y
);
hold on
;
X
=linspace
(3,4,100
);
Y
=f
(3
)+
(f
(4
)-f
(3
))/
(x
(4
)-x
(3
)).*
(X-x
(3
))-
(M
(4
)/6+M
(3
)/3
)*
(x
(4
)-x
(3
)).*
(X-x
(3
))+M
(3
)/2.*
(X-x
(3
)).^2+
(M
(4
)-M
(3
))/
(x
(4
)-x
(3
)).*
(X-x
(3
)).^3/6
;
plot
(X,Y
);
hold on
;
X
=linspace
(4,6,200
);
Y
=f
(4
)+
(f
(5
)-f
(4
))/
(x
(5
)-x
(4
)).*
(X-x
(4
))-
(M
(5
)/6+M
(4
)/3
)*
(x
(5
)-x
(4
)).*
(X-x
(4
))+M
(4
)/2.*
(X-x
(4
)).^2+
(M
(5
)-M
(4
))/
(x
(5
)-x
(4
)).*
(X-x
(4
)).^3/6
;
plot
(X,Y
);
hold on
;
legend
('[1,2]',
'[2,3]',
'[3,4]',
'[4,6]');