在许多问题中,通常根据实验、观测或经验得到的函数表或离散点上的信息,去研究分析函数的有关特性。其中插值法是一种最基本的方法,以下给出最基本的插值问题一三次样条插值的基本方法:
对插值区间[a,b]进行划分,a<x0<x1<…<xn≤b.函数y=f(x)在节点xi上的值yi=f(xi)(i=0.,1.2…n),并且如果函数S(x)在每个小区间[xi,x(i+1)]上是三次多项式,在[a,b]上有二阶连续导数,则称S(x)是[a,b]上的三次样条函数,如果S(x)在节点xi上还满足条件S(xi)=yi(i=0,1…n),则称S(x)为三次样条插值函数。
解 :根据题意,编写matlab代码如下:
clear all
clc
x = [0.25 0.3 0.39 0.45 0.53]
y = [0.5000 0.5477 0.6245 0.6708 0.7280]
n= length(x);
for i=1:n- 1
h(i)=x(i+ 1)-x(i);
end
for i= 1:n- 2
k(i)=h(i+1)/(h(i) +h(i+ 1));
u(i)= h(i)/(h(i) + h(i+ 1));
end
for i = 1:n-2
g1(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1) +k(i)*(y(i+ 1)- y(i))/h(i));
end
g0=3*(y(2)- y(1))/h(1);
g00=3* (y(n)- y(n- 1))/h(n- 1);
g=[g0 g1 g00];
g= transpose(g)
k1 = [k 1];
u1 = [1 u];
Q = 2*eye(5) + diag(u1,1) + diag(k1,-1)
m = transpose(Q\g)
syms X;
for i=1:n-1
p1(i)=(1+2* (X-x(i))/h(i))*((X-x(i+ 1))/h(i))^2* y(i);
p2(i)=(1-2* (X-x(i+1))/h(i))* ((X- x(i))/h(i))^2*y(i+ 1);
p3(i)= (X-x(i)) * ((X-x(i+ 1))/h(i))^2* m(i);
p4(i)= (X-x(i+ 1)) * ((X- x(i))/h(i))^2*m(i+ 1);
p(i)=p1(1)+p2(1) + p3(1) + p4(1);
end
s1 = p(1)
s2 = p(2)
s3 = p(3)
s4 = p(4)
for k=1:4
for z = x(k):0.001:x(k+ 1)
q = eval(subs(p(k), 'X','z'));
end
end
运行结果如下:
这里的结果可以自己运行一下。