数值分析三次样条插值确定边界条件的函数表达式求解(MATLAB)实验报告(附实验代码)

一.实验内容
实现三次样条插值,给定从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)=2
fn__;
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)/(6
h(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(这是第二个问题的数据)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值