目录
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,可以忽略不计。