系列1 详解matlab实现三次样条插值的三弯矩法

目录

1 习题一

1.1 题目

1.2 问题背景

1.3 算法

1.4 matlab编程实现

1.5 结果展示

1.6 讨论

---

1 习题一

1.1 题目

已知如下数据：

x	-1.00	-0.54	0.13	1.12	1.89	2.06	2.54	2.82	3.50
f(x)	-2.46	-5.26	-1.87	0.05	1.65	2.69	4.56	7.89	10.31

且端点约束条件为𝑓^′ (−1)=5,𝑓^′ (3.50)=29.16，

用三次样条插值的三弯矩法求𝑓(−0.02)和𝑓(2.56)。

1.2 问题背景

分段低次插值可以说很好地克服了Runge 现象，并且表达式也很简单，但是它们的光滑性不够理想，分段线性插值只具有连续性，分段三次Hermite插值才是一次连续可微的。这往往不能满足实际问题所提出的光滑性要求。

在工业设计中，对曲线的光滑性均有一定的要求,例如，在飞机、船舶、汽车等外形设计中，要求外形曲线呈流线型，早期为设计这种具有流线型的曲线，通常采用样条这种绘图工具,它是一根富有弹性的细长木条,工程师把它用压铁固定在已知点处转折自如,这种曲线称为样条曲线。其中最常用的是三次样条曲线,三弯矩法是力学问题中常用的三次样条曲线解法。

1.3 算法

1.4 matlab编程实现

（将函数文件放在与主程序同一个文件夹下，然后运行主程序）

主程序：

clear;

x=[-1.00  -0.54  0.13  1.12  1.89  2.06  2.54  2.82  3.50];

y=[-2.46  -5.26 -1.87  0.05  1.65  2.69  4.56  7.89  10.31];

f0=5;fn=29.16;

x0=[-0.02 ,2.56];

  % x为插值节点，y为插值节点处的函数值，f0为f'（x(1)），fn为f'（x(n)）,x0为待求插值点或矩阵

y3=sanwanju(x,y,f0,fn,x0)  %统一调用，y3(1)≈f(-0.02);y3(2)≈f(2.56)

% y1=sanwanju(x,y,f0,fn,-0.02)  %等价形式：分别调用求f(-0.02)的近似值

% y2=sanwanju(x,y,f0,fn,2.56)   %分别调用求f(2.56)的近似值

函数1(三弯矩法求指定点三次样条插值近似值)：

function y0=sanwanju(x,y,f0,fn,x0)

%用第一种边界条件的三弯矩方程组,来求三次样条插值函数

% x为插值节点，y为插值节点处的函数值，f0为f'（x(1)），fn为f'（x(n)）,x0为待求插值点，y0为插值近似值

syms w;

n=length(x);

a=length(x0);

h=zeros(n-1,1);f=zeros(n-1,1);d=zeros(n,1);

for i=1:n-1

    h(i)=x(i+1)-x(i);  %存储步长

    f(i)=(y(i+1)-y(i))/h(i);  %存储二阶差商

end

d(1)=6/h(1)*(f(1)-f0);

d(n)=6/h(n-1)*(fn-f(n-1));

u=zeros(n-2,1);lamda=zeros(n-2,1);

for i=1:n-2

    u(i)=h(i)/(h(i)+h(i+1));

    lamda(i)=1-u(i);

    d(i+1)=6/(h(i)+h(i+1))*(f(i+1)-f(i));

end

A=2*eye(n,n);A(1,2)=1;A(n,n-1)=1;  %将计算所得赋值给A矩阵对应位置

for i=1:n-2

    A(i+1,i)=u(i);

    A(i+1,i+2)=lamda(i);

end

M=chase(A,d);  %用追赶法求解M

for i=1:a

    for j=1:n-1

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

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

            S=@(w)M(j)*((x(j+1)-w)^3/(6*h))+M(j+1)*((w-x(j))^3/(6*h))+(y(j)-M(j)*h*h/6)*(x(j+1)-w)/h+(y(j+1)-M(j+1)*h*h/6)*(w-x(j))/h;        

            y0(i)=feval(S,x0(i));  %相当于输出y0(i)=S（x0(i)）

            break;

        end

    end

end

函数2（用追赶法求解三弯矩方程组）：

function x=chase(A,b)

%追赶法求解线性方程组

n=length(A);

L=eye(n,n);U=eye(n,n);z=zeros(n,1);x=zeros(n,1);

for k=1:n-1

    U(k,k+1)=A(k,k+1);  %赋给U矩阵元素c(i)值，与A矩阵相同

end

for i=2:n

    U(1,1)=A(1,1);  %U矩阵赋初值

    L(i,i-1)=A(i,i-1)/U(i-1,i-1);  %逐个计算L矩阵元素l(i)值,其中l(i)=L(i,i-1)

    U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);  %逐个计算U矩阵元素r(i)值,其中r(i)=U(i,i)

end

% z=inv(L)*b;

% x=inv(U)*z;

z(1)=b(1); 

for i=2:n

    %回代求解z=inv(L)*b

    z(i)=b(i)-L(i,i-1)*z(i-1);

end

x(n)=z(n)/U(n,n);   

for i=n-1:-1:1

    %回代求解x=inv(U)*z

    sum=0;

    for j=i+1:n

        sum=sum+U(i,j)*x(j);

    end

    x(i)=(z(i)-sum)/U(i,i);

end

end

1.5 结果展示

>> work1

yy =

   -3.1058    4.7834

1.6 讨论

csape与spline函数比较相近，都可以选择样条的边界条件，interp1无法使用边界条件；而相较于csape，spline函数，interp1可以选择几种不同的插值方法。由于本题是一阶导数边界条件，下面分别选择csape与spline函数进行验证。

主程序添加下列程序：

S=spline(x,[f0,y,fn]);yy2=ppval(S,x0)  %spline函数验证

error1=yy-yy2

cs= csape(x,[f0,y,fn],[1,1]);yy3=ppval(cs,x0)  %csape函数验证

error2=yy-yy3

再次运行：

>> work1

yy =

   -3.1058    4.7834

yy2 =

   -3.1058    4.7834

error1 =

   1.0e-15 *

    0.4441    0.8882

yy3 =

   -3.1058    4.7834

error2 =

   1.0e-14 *

    0.1332    0.0888

可以看到error1与error2极小，它们的误差量级分别为1.0e-15，1.0e-14，可以忽略不计。