数值分析第一次作业.docx
《数值分析第一次作业.docx》由会员分享，可在线阅读，更多相关《数值分析第一次作业.docx（26页珍藏版）》请在冰点文库上搜索。
数值分析第一次作业
问题1：
20.给定数据如下表：
xj
0.25
0.30
0.39
0.45
0.53
yj
0.5000
0.5477
0.6245
0.6708
0.7280
试求三次样条插值S（x），并满足条件
（1）S`（0.25）=1.0000，S`（0.53）=0.6868；
（2）S’’（0.25）=S’’（0.53）=0。
分析：
本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，
（1）是已知一阶倒数，
（2）是已知自然边界条件。
对于第一种边界（已知边界的一阶倒数值），可写出下面的矩阵方程。
其中
j=
，
i=
，dj=6f[xj-1,xj,xj+1]，
n=1，
0=1
对于第一种边界条件d0=
（f[x0,x1]-f0`），dn=
（f`n-f`[xn-1,xn]）
解：
由matlab计算得：
x
y
h
d
0.25
0.5000
-5.5200
0.30
0.5477
0.0500
0.3571
1
-4.3143
0.39
0.6245
0.0900
0.6000
0.6429
-3.2667
0.45
0.6708
0.0600
0.4286
0.4000
-2.4286
0.53
0.7280
0.0800
1.000
0.5714
-2.1150
由此得矩阵形式的线性方程组为：
解得M0=-2.0286;M1=-1.4627;M2=-1.0333;M3=-0.8058;M4=-0.6546
S（x）=
Matlab程序代码如下：
functiontgsanci（n,s,t）%n代表元素数，s,t代表端点的一阶导。
x=[0.250.300.390.450.53];
y=[0.50.54770.62450.67080.7280];
n=5,s=1.0,t=0.6868
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=6*（f
（1）-s）/h
（1）
d（n）=6*（t-f（n-1））/h（n-1）
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=1;
u（n）=1;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
m=b*d'
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：
其中
j=
，
i=
，dj=6f[xj-1,xj,xj+1]，
n=
0=0d0=dn=0
解：
由matlab计算得：
x
y
h
dn
0.25
0.5000
0
0.30
0.5477
0.05
0.3571
0
-4.3143
0.39
0.6245
0.09
0.6000
0.6429
-3.2667
0.45
0.6708
0.06
0.4286
0.4000
-2.4268
0.53
0.7280
0.08
0
0.5714
0
由此得矩阵形式的线性方程组为：
解得M0=0;M1=-1.8795;M2=-0.8636;M3=-1.0292;M4=0
S（x）=
matlab程序代码如下：
functiontgsanci（n,s,t）%n代表元素数，
x=[0.250.300.390.450.53];
y=[0.50.54770.62450.67080.7280];
n=5
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=0
d（n）=0
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=0;
u（n）=0;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
k=a
m=b*d'
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
p
由下面的一段程序，画出自然边界与第一边界的图形：
程序代码如下：
x=[0.250.300.390.450.53];
y=[0.50.54770.62450.67080.7280];
s=csape（x,y,'variational'）
fnplt（s,'r'）
holdon
xlabel（'x'）
ylabel（'y'）
gtext（'三次样条自然边界'）
plot（x,y,'o',x,y,':
g'）
holdon
s.coefs
r=csape（x,y,'complete',[1.00.6868]）
fnplt（r,'b'）
gtext（'三次样条第一边界'）
holdon
r.coefs
图中的三条线几乎重合。
图20-1在一小段区间内的图形
图20-2在整个给出的区间的图形
问题2：
1已知函数在下列各点的值为
xi
0.2
0.4
0.6
.0.8
1.0
f（xi）
0.98
0.92
0.81
0.64
0.38
试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。
用图给出{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}，P4（x）及S（x）。
分析：
（1）要用4次牛顿插值多项式处理数据，
Pn=f（x0）+f[x0,x1]（x-x0）+f[x0,x1,x2]（x-x0）（x-x1）+···+f[x0,x1，···xn]（x-x0）···（x-xn-1）
首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。
解：
由matlab计算得：
xi
f（xi）
一阶差商
二阶差商
三阶差商
四阶差商
0.20
0.980
0.40
0.920
-0.30000
0.60
0.810
-0.55000
-0.62500
0.80
0.640
-0.85000
-0.75000
-0.20833
1.00
0.380
-1.30000
-1.12500
-0.62500
-0.52083
所以有四次插值牛顿多项式为：
P4（x）=0.98-0.3（x-0.2）-0.62500（x-0.2）（x-0.4）-0.20833（x-0.2）（x-0.4）（x-0.6）-0.52083（x-0.2）（x-0.4）（x-0.6）（x-0.8）
计算牛顿插值中多项式系数的程序如下：
functionvarargout=newtonliu（varargin）
clear,clc
x=[0.20.40.60.81.0];
fx=[0.980.920.810.640.38];
newtonchzh（x,fx）;
functionnewtonchzh（x,fx）
%由此函数可得差分表
n=length（x）;
fprintf（'*****************差分表*****************************\n'）;
FF=ones（n,n）;
FF（:
1）=fx';
fori=2:
n
forj=i:
n
FF（j,i）=（FF（j,i-1）-FF（j-1,i-1））/（x（j）-x（j-i+1））;
end
end
fori=1:
n
fprintf（'%4.2f',x（i））;
forj=1:
i
fprintf（'%10.5f',FF（i,j））;
end
fprintf（'\n'）;
end
（2）用三次样条插值函数由上题分析知,要求各点的M值：
有matlab编程计算求出个三次样条函数：
解得：
M0=0;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=0
三次样条插值函数计算的程序如下：
functiontgsanci（n,s,t）%n代表元素数，s,t代表端点的一阶导。
x=[0.20.40.60.81.0];
y=[0.980.920.810.640.38];
n=5
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=0
d（n）=0
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=0;
u（n）=0;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
m=b*d'
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
p
下面是画牛顿插值以及三次样条插值图形的程序：
x=[0.20.40.60.81.0];
y=[0.980.920.810.640.38];
plot（x,y）
holdon
fori=1:
1:
5
y（i）=0.98-0.3*（x（i）-0.2）-0.62500*（x（i）-0.2）*（x（i）-0.4）-0.20833*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）-0.52083*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）*（x（i）-0.8）
end
k=[011011]
x0=0.2+0.08*k
fori=1:
1:
4
y0（i）=0.98-0.3*（x（i）-0.2）-0.62500*（x（i）-0.2）*（x（i）-0.4）-0.20833*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）-0.52083*（x（i）-0.2）*（x（i）-0.4）*（x（i）-0.6）*（x（i）-0.8）
end
plot（x0,y0,'o',x0,y0）
holdon
y1=spline（x,y,x0）
plot（x0,y1,'o'）
holdon
s=csape（x,y,'variational'）
fnplt（s,'r'）
holdon
gtext（'三次样条自然边界'）
gtext（'原图像'）
gtext（'4次牛顿插值'）
运行上述程序可知：
给出的{（xi，yi），xi=0.2+0.08i，i=0，1,11,10}点，S（x）及P4（x）图形如下所示：
问题3：
3下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9各点作8次多项式插值L8（x）.
（2）用三次样条（自然边界条件）程序求S（x）。
从结果看在[0,64]上，那个插值更精确；在区间[0,1]上，两种哪个更精确？
分析：
L8（x）可由公式Ln（x）=
得出。
三次样条可以利用自然边界条件。
写成矩阵：
其中
j=
，
i=
，dj=6f[xj-1,xj,xj+1]，
n=
0=0d0=dn=0
解：
l0（x）=
l1（x）=
l2（x）=
l3（x）=
l4（x）=
l5（x）=
l6（x）=
l7（x）=
l8（x）=
L8（x）=l1（x）+2l2（x）+3l3（x）+4l4（x）+5l5（x）+6l6（x）+7l7（x）+8l8（x）
求三次样条插值函数由matlab计算：
可得矩阵形式的线性方程组为：
解得:
M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=--0.0009;
M8=0
S（x）=
拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。
图3-1为[01]的曲线，图3-2为[064]区间上的曲线。
由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在[01]朗格朗日插值更精确。
而在区间[064]上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在[3070]之间有很大的振荡，所在在区间[064]三次样条插值更精确写。
图3-1
图3-2
Matlab程序代码如下：
求三次样条插值函数的程序：
functiontgsanci（n,s,t）%n代表元素数，s,t代表端点的一阶导。
y=[012345678];
x=[01491625364964];
n=9
forj=1:
1:
n-1
h（j）=x（j+1）-x（j）;
end
forj=2:
1:
n-1
r（j）=h（j）/（h（j）+h（j-1））;
end
forj=1:
1:
n-1
u（j）=1-r（j）;
end
forj=1:
1:
n-1
f（j）=（y（j+1）-y（j））/h（j）;
end
forj=2:
1:
n-1
d（j）=6*（f（j）-f（j-1））/（h（j-1）+h（j））;
end
d
（1）=0
d（n）=0
a=zeros（n,n）;
forj=1:
1:
n
a（j,j）=2;
end
r
（1）=0;
u（n）=0;
forj=1:
1:
n-1
a（j+1,j）=u（j+1）;
a（j,j+1）=r（j）;
end
b=inv（a）
m=b*d'
t=a
p=zeros（n-1,4）;%p矩阵为S（x）函数的系数矩阵
forj=1:
1:
n-1
p（j,1）=m（j）/（6*h（j））;
p（j,2）=m（j+1）/（6*h（j））;
p（j,3）=（y（j）-m（j）*（h（j）^2/6））/h（j）;
p（j,4）=（y（j+1）-m（j+1）*（h（j）^2/6））/h（j）;
end
p
%画图形比较那个插值更精确的函数：
x0=[01491625364964];
y0=[012345678];
x=0:
64;y=lagr1（x0,y0,x）;
plot（x0,y0,'o'）
holdon
plot（x,y,'r'）;
holdon;
pp=csape（x0,y0,'variational'）
fnplt（pp,'g'）;
holdon;
plot（x0,y0,':
b'）;holdon
%axis（[0201]）;%看[01]区间的图形时加上这条指令
gtext（'三次样条插值'）
gtext（'原图像'）
gtext（'拉格朗日插值'）
%下面是求拉格朗日插值的函数
functiony=lagr1（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
问题：
16.观测物体的直线运动，得出以下数据：
时间（t/s）
0
0.9
1.9
3.0
3.9
5.0
距离（s/m）
0
10
30
50
80
110
求运动方程。
分析：
由所给的数据在坐标纸上描出如图16-1所示。
从图中可以看到各点在一条直线的附近，故可以选择线性函数作拟合曲线，即令s（t）=a+bt,这里m=5,n=1。
法方程矩阵为下面的形式：
的线性无关性，知道该方程存在唯一解
解：
由上面的分析以及通过matlab计算得：
法方程矩阵如下：
解之得：
a=-7.8550,b=22.2538。
由此可得运动方程为：
S（t）=22.2538t-7.8550.
在matlab中拟合的曲线如图16-2所示：
图16-1图16-2
Matlab源程序代码：
计算部分的程序代码：
x=[00.91.93.03.95.0];
y=[010305080110];
r=zeros（2,2）;
fori=1:
1:
6;
r（1,1）=r（1,1）+1;
end
fori=1:
1:
6
r（1,2）=r（1,2）+x（i）;
end
fori=1:
1:
6
r（2,1）=r（2,1）+x（i）;
end
fori=1:
1:
6
r（2,2）=r（2,2）+x（i）^2;
end
a=zeros（2,1）;
fori=1:
1:
6
a（1,1）=a（1,1）+y（i）;
a（2,1）=a（2,1）+y（i）*x（i）
end
k=r
t=a
d=inv（r）;a=d*a
画图的程序代码：
x=[00.91.93.03.95.0];
y=[010305080110];
xx=0:
0.02:
5.0;
pp=polyfit（x,y,1）;
yy=polyval（pp,xx）;
plot（x,y,'o',xx,yy）
问题：
18在某化学反应中，由试验得分解物浓度与时间的关系如下：
时间（t/s）
0
5
10
15
20
25
30
35
40
45
50
55
浓度y/（10-4）
0
1.27
2.16
2.86
3.44
3.87
4.15
4.37
4.51
4.58
4.62
4.64
用最小二乘法求y=f（t）。
分析：
首先要选择拟合曲线，作出给定数据的散点图如下：
图18-1给定数据的散点图
通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。
具有这种特性的曲线很多，我们选取如下三种函数来拟合：
①多项式（x）=a0a1x…amxm，m为适当选取的正整数；
②
；
③
（a>0,b<0）。
在matlab中分别用上述拟合函数拟合曲线，本题应采用倒指数拟合，即y（t）=a*exp（bt-1）拟合曲线。
对y（t）=a*exp（bt-1）两边取对数有㏑y=㏑a+b/t;如令Y=
㏑y，A=㏑a，则得Y=A+b/t;为确定A，b,先将（ti，yi）转化为（ti,Yi）.
根据最小二乘法，取
，
。
用lsqcurvefit函数拟合曲线，程序代码如下：
创建M文件：
functionF=mf（a,t）
F=a
（1）*exp（a
（2）.*t.^（-1））;
在命令窗口中敲入如下代码：
t=[5:
5:
55]
y=[1.272.162.863.443.874.154.374.514.584.624.64];
a0=[0.5,0.5];
a=lsqcurvefit（'pf',a0,t,y）
回车后可得：
a
（1）=5.5031,a
（2）=-8.7872
下面的计算问题用matlab编程实现：
x=[510152025303540455055];
y0=[1.272.162.863.443.874.154.374.514.584.624.64];
y=log（y0）
y
（1）=0
r=zeros（2,2）;
fori=1:
1:
12;
r（1,1）=r（1,1）+1;
end
fori=2:
1:
12
r（1,2）=r（1,2）+1/x（i）;
end
fori=2:
1:
12
r（2,1）=r（2,1）+1/x（i）;
end
fori=2:
1:
12
r（2,2）=r（2,2）+1/x（i）^2;
end
a=zeros（2,1）;
f
升级会员 