一.实验内容
实现三次样条插值,给定从Xo到Xn的点,再给定边界条件,运用数学方法求出该三次样条函数。其中,边界条件有两种,第一种:给定边界的一阶导数;第二种:给定边界的二阶导数。
二:实验工具
MATLAB
三.实验思路
实验的开头,用load()函数输入数据点,这样做的目的是使输入数据方便快捷,load()函数从MATLAB文件所在地方读取data.txt文件,分别将x和y存在两个列向量中,其中,当我们往data.txt文件中输入数据时,注意X在前Y在后,且二者中间隔一个空格,另外,两个相邻的数据之间用回车键分割,还要注意数据对齐,以免引起不必要的麻烦。
实验多采用for循环来分别计算hi, di, f[xi,xi+1](代码中以g(i)表示)等等然后分别将各数存在列向量中,以便后面使用
采用对话框的方式,输出用户可以选择的作为输入的边界条件,然后用switch语句,分别操作。
具体算函数的形式还有矩阵的形式时,参考课本的推到过程:
四.实验代码
a=load(“data.txt”)
x = a(:,1);
y = a(:,2);
n=input(‘输入n:’);
A=zeros(n+1,n+1);
for i=1:n+1
A(i,i)=2;
end%定义矩阵A对角线上的元素
for i=1:n
h(i)=x(i+1)-x(i);
end%定义hi,因MATLAB中不存在以0为下标,因此存储1到n的h
for i=2:n
v(i)=h(i)/(h(i-1)+h(i));
end
for i=2:n
u(i)=h(i-1)/(h(i-1)+h(i));
end
for i=1:n
g(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
end%定义各一阶差商
for i=2:n
d(i)=(6*(g(i)-g(i-1)))/(h(i-1)+h(i));
end%定义di
for i=2:n
k(i)=d(i);
end
for i=2:n
A(i,i+1)=v(i);
A(i,i-1)=u(i);
end
fprintf(‘1:已知X0和Xn处的一阶导数\n2:已知X0和Xn处的二阶倒数\n’)
sel=input(‘根据以上信息输入边界条件类型’);
switch sel
case 2
f0__=input('输入X0处的二阶倒数 ‘);
fn__=input(‘输入Xn处的二阶倒数’);
v(1)=0;
A(1,2)=v(1);
u(n+1)=0;
A(n+1,n)=u(n+1);
k(1)=2f0__;
k(n+1)=2fn__;
case 1
f0_=input(‘ 输入X0处的一阶倒数’);
fn_=input(‘输入Xn处的一阶倒数’);
v(1)=1;
A(1,2)=v(1);
k(1)=(6*(g(1)-f0_))/h(1);
u(n+1)=1;
A(n+1,n)=u(n+1);
k(n+1)=(6*(fn_-g(n)))/h(n);
end
for i=1:n+1
B(i,1)=k(i)
end
fprintf(‘解得M分别为:’)
M=A\B
for i=1:n
fprintf('当%f≤x≤%f时,函数的表达式为: ',x(1),x(2));
syms X
syms Y
X=x(i):0.01:x(i+1);
a1=M(i)/(6h(i));
a2=M(i+1)/(6h(i));
b1=(y(i)-M(i)*h(i)h(i)/6)/h(i);
b2=(y(i+1)-M(i+1)h(i)h(i)/6)/h(i);
fprintf(‘S(x)=%f(%f-X)3+%f(X-%f)3+%f(%f-X)+%f(X-%f)’,a1,x(i+1),a2,x(i),b1,x(i+1),x(i));
Y=a1.(x(i+1)-X).3+a2.*(X-x(i)).3+b1.(x(i+1)-X)+b2.(X-x(i))
plot(X,Y)
hold on
end
hold off
五.实验结果:
①以课本课后习题20为例:
第一问:
a =
0.2500 0.5000
0.3000 0.5477
0.3900 0.6245
0.4500 0.6708
0.5300 0.7280
输入n:4
1:已知X0和Xn处的一阶倒数
2:已知X0和Xn处的二阶倒数
根据以上信息输入边界条件类型:1
输入X0处的一阶倒数1
输入Xn处的一阶倒数0.6868
解得M分别为:
M =
-2.0286
-1.4627
-1.0333
-0.8058
-0.6546
当0.250000≤x≤0.300000时,函数的表达式为: S(x)=-6.762098(0.300000-X)3±4.875803(X-0.250000)3+10.016905(0.300000-X)+0.250000(X-
Y =
0.5000 0.5099 0.5196 0.5291 0.5385 0.5477
当0.300000≤x≤0.390000时,函数的表达式为: S(x)=-2.708780(0.390000-X)3±1.913602(X-0.300000)3+6.107497(0.390000-X)+0.300000(X-
Y =
0.5477 0.5568 0.5657 0.5744 0.5831 0.5916 0.6000 0.6083 0.6164 0.6245
当0.390000≤x≤0.450000时,函数的表达式为: S(x)=-2.870403(0.450000-X)3±2.238418(X-0.390000)3+10.418667(0.450000-X)+0.390000(X-
Y =
0.6245 0.6325 0.6403 0.6481 0.6557 0.6633 0.6708
当0.450000≤x≤0.530000时,函数的表达式为: S(x)=-1.678813(0.530000-X)3±1.363718(X-0.450000)3+8.395744(0.530000-X)+0.450000(X-
Y =
0.6708 0.6782 0.6855 0.6928 0.7000 0.7071 0.7141 0.7211 0.7280
第二问:
a =
0.2500 0.5000
0.3000 0.5477
0.3900 0.6245
0.4500 0.6708
0.5300 0.7280
输入n:4
1:已知X0和Xn处的一阶倒数
2:已知X0和Xn处的二阶倒数
3:f(x)是以Xn-X0为周期的周期函数
根据以上信息输入边界条件类型:2
输入X0处的二阶倒数0
输入Xn处的二阶倒数0
解得M分别为:
M =
0
-1.8795
-0.8636
-1.0292
0
当0.250000≤x≤0.300000时,函数的表达式为: S(x)=0.000000(0.300000-X)3±6.265165(X-0.250000)3+10.000000(0.300000-X)+0.250000(X-
Y =
0.5000 0.5097 0.5193 0.5289 0.5384 0.5477
当0.300000≤x≤0.390000时,函数的表达式为: S(x)=-3.480647(0.390000-X)3±1.599303(X-0.300000)3+6.113749(0.390000-X)+0.300000(X-
Y =
0.5477 0.5568 0.5658 0.5746 0.5832 0.5917 0.6001 0.6083 0.6165 0.6245
当0.390000≤x≤0.450000时,函数的表达式为: S(x)=-2.398955(0.450000-X)3±2.858954(X-0.390000)3+10.416970(0.450000-X)+0.390000(X-
Y =
0.6245 0.6324 0.6403 0.6481 0.6557 0.6633 0.6708
当0.450000≤x≤0.530000时,函数的表达式为: S(x)=-2.144216(0.530000-X)3+0.000000(X-0.450000)3+8.398723(0.530000-X)+0.450000(X-
Y =
0.6708 0.6782 0.6855 0.6927 0.6998 0.7069 0.7140 0.7210 0.7280
②以课本44页例7为例:
test3
a =
27.7000 4.1000
28.0000 4.3000
29.0000 4.1000
30.0000 3.0000
输入n:3
1:已知X0和Xn处的一阶倒数
2:已知X0和Xn处的二阶倒数
3:f(x)是以Xn-X0为周期的周期函数
根据以上信息输入边界条件类型:1
输入X0处的一阶倒数3
输入Xn处的一阶倒数-4
解得M分别为:
M =
-23.5314
0.3960
0.8297
-9.1149
当27.700000≤x≤28.000000时,函数的表达式为: S(x)=-13.072974(28.000000-X)3+0.220022(X-27.700000)3+14.843234(28.000000-X)+27.700000(X-
Y =
列 1 至 13
4.1000 4.1288 4.1554 4.1798 4.2020 4.2222 4.2405 4.2569 4.2715 4.2844 4.2956 4.3053 4.3135
列 14 至 26
4.3204 4.3259 4.3301 4.3332 4.3353 4.3363 4.3364 4.3357 4.3342 4.3321 4.3293 4.3261 4.3223
列 27 至 31
4.3183 4.3139 4.3094 4.3047 4.3000
当28.000000≤x≤29.000000时,函数的表达式为: S(x)=0.066007(29.000000-X)3+0.138284(X-28.000000)3+4.233993(29.000000-X)+28.000000(X-
Y =
列 1 至 13
4.3000 4.2953 4.2907 4.2861 4.2815 4.2770 4.2725 4.2681 4.2637 4.2593 4.2550 4.2508 4.2465
列 14 至 26
4.2424 4.2382 4.2342 4.2301 4.2261 4.2222 4.2183 4.2144 4.2106 4.2069 4.2032 4.1995 4.1959
列 27 至 39
4.1924 4.1889 4.1854 4.1820 4.1787 4.1754 4.1722 4.1690 4.1658 4.1628 4.1597 4.1568 4.1538
列 40 至 52
4.1510 4.1482 4.1454 4.1428 4.1401 4.1376 4.1351 4.1326 4.1302 4.1279 4.1256 4.1234 4.1212
列 53 至 65
4.1192 4.1171 4.1152 4.1133 4.1114 4.1097 4.1079 4.1063 4.1047 4.1032 4.1018 4.1004 4.0991
列 66 至 78
4.0978 4.0966 4.0955 4.0945 4.0935 4.0926 4.0918 4.0910 4.0903 4.0897 4.0892 4.0887 4.0883
列 79 至 91
4.0879 4.0877 4.0875 4.0874 4.0874 4.0874 4.0875 4.0877 4.0880 4.0883 4.0887 4.0892 4.0898
列 92 至 101
4.0905 4.0912 4.0920 4.0929 4.0939 4.0950 4.0961 4.0973 4.0986 4.1000
当29.000000≤x≤30.000000时,函数的表达式为: S(x)=0.138284(30.000000-X)3±1.519142(X-29.000000)3+3.961716(30.000000-X)+29.000000(X-
Y =
列 1 至 13
4.1000 4.1015 4.1030 4.1046 4.1063 4.1080 4.1097 4.1114 4.1132 4.1150 4.1167 4.1185 4.1202
列 14 至 26
4.1219 4.1235 4.1251 4.1266 4.1281 4.1294 4.1307 4.1318 4.1329 4.1338 4.1346 4.1352 4.1357
列 27 至 39
4.1360 4.1361 4.1361 4.1358 4.1354 4.1347 4.1338 4.1327 4.1313 4.1297 4.1278 4.1256 4.1231
列 40 至 52
4.1204 4.1173 4.1140 4.1103 4.1062 4.1019 4.0971 4.0920 4.0866 4.0807 4.0745 4.0678 4.0608
列 53 至 65
4.0533 4.0453 4.0370 4.0282 4.0189 4.0091 3.9989 3.9881 3.9769 3.9651 3.9529 3.9400 3.9267
列 66 至 78
3.9128 3.8983 3.8833 3.8676 3.8514 3.8346 3.8171 3.7991 3.7804 3.7610 3.7411 3.7204 3.6991
列 79 至 91
3.6771 3.6544 3.6310 3.6068 3.5820 3.5564 3.5301 3.5031 3.4752 3.4466 3.4172 3.3871 3.3561
列 92 至 101
3.3243 3.2917 3.2582 3.2240 3.1888 3.1528 3.1159 3.0782 3.0395 3.0000
我的data.txt的内容是:
27.7 4.1
28.0 4.3
29.0 4.1
30.0 3.0(这是第二个问题的数据)