数值分析第一次作业.doc

内科大数值分析作业
班级：研5班 姓名：王雪萍 学号：201102208 专业：双控
-----------------------------------------------------------------------------------------
问题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程序代码如下：
function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
x=[0.25 0.30 0.39 0.45 0.53];
y=[0.5 0.5477 0.6245 0.6708 0.7280];
n=5,s=1.0,t=0.6868
for j=1:1:n-1
h(j)=x(j+1)-x(j);
end
for j=2:1:n-1
r(j)=h(j)/(h(j)+h(j-1));
end
for j=1:1:n-1
u(j)=1-r(j);
end
for j=1:1:n-1
f(j)=(y(j+1)-y(j))/h(j);
end
for j=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);
for j=1:1:n
a(j,j)=2;
end
r(1)=1;
u(n)=1;
for j=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)函数的系数矩阵
for j=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=0 d0=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程序代码如下：
function tgsanci(n,s,t) %n代表元素数，
x=[0.25 0.30 0.39 0.45 0.53];
y=[0.5 0.5477 0.6245 0.6708 0.7280];
n=5
for j=1:1:n-1
h(j)=x(j+1)-x(j);
end
for j=2:1:n-1
r(j)=h(j)/(h(j)+h(j-1));
end
for j=1:1:n-1
u(j)=1-r(j);
end
for j=1:1:n-1
f(j)=(y(j+1)-y(j))/h(j);
end
for j=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);
for j=1:1:n
a(j,j)=2;
end
r(1)=0;
u(n)=0;
for j=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)函数的系数矩阵
for j=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.25 0.30 0.39 0.45 0.53];
y=[0.5 0.5477 0.6245 0.6708 0.7280];
s=csape(x,y,'variational')
fnplt(s,'r')
hold on
xlabel('x')
ylabel('y')
gtext('三次样条自然边界')
plot(x,y,'o',x,y,':g')
hold on
s.coefs
r=csape(x,y,'complete',[1.0 0.6868])
fnplt(r,'b')
gtext('三次样条第一边界')
hold on
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）
计算牛顿插值中多项式系数的程序如下：
function varargout=newtonliu(varargin)
clear,clc
x=[0.2 0.4 0.6 0.8 1.0];
fx=[0.98 0.92 0.81 0.64 0.38];
newtonchzh(x,fx);
function newtonchzh(x,fx)
%由此函数可得差分表
n=length(x);
fprintf('*****************差分表*****************************\n');
FF=ones(n,n);
FF(:,1)=fx';
for i=2:n
for j=i:n
FF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));
end
end
for i=1:n
fprintf('%4.2f',x(i));
for j=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
三次样条插值函数计算的程序如下：
function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
x=[0.2 0.4 0.6 0.8 1.0];
y=[0.98 0.92 0.81 0.64 0.38];
n=5
for j=1:1:n-1
h(j)=x(j+1)-x(j);
end
for j=2:1:n-1
r(j)=h(j)/(h(j)+h(j-1));
end
for j=1:1:n-1
u(j)=1-r(j);
end
for j=1:1:n-1
f(j)=(y(j+1)-y(j))/h(j);
end
for j=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);
for j=1:1:n
a(j,j)=2;
end
r(1)=0;
u(n)=0;
for j=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)函数的系数矩阵
for j=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.2 0.4 0.6 0.8 1.0];
y=[0.98 0.92 0.81 0.64 0.38];
plot(x,y)
hold on
for i=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=[0 1 10 11]
x0=0.2+0.08*k
for i=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 )
hold on
y1=spline(x,y,x0)
plot(x0,y1,'o')
hold on
s=csape(x,y,'variational')
fnplt(s,'r')
hold on
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=0 d0=dn=0
解：l0(x)=
l1(x)=
l2(x)=
l3(x)=
l4(x)=
l5(x)=
l6(x)=
l7(x)=
l8(x)=
L8(x)= l1(x)+2 l2(x)+3 l3(x)+4 l4(x)+5 l5(x)+6 l6(x)+7 l7(x)+8 l8(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为[0 1]的曲线，图3-2为[0 64]区间上的曲线。由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在[0 1]朗格朗日插值更精确。而在区间[0 64]上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在[30 70]之间有很大的振荡，所在在区间[0 64]三次样条插值更精确写。

图3-1

图3-2
Matlab程序代码如下：
求三次样条插值函数的程序：
function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
y=[0 1 2 3 4 5 6 7 8];
x=[0 1 4 9 16 25 36 49 64];
n=9
for j=1:1:n-1
h(j)=x(j+1)-x(j);
end
for j=2:1:n-1
r(j)=h(j)/(h(j)+h(j-1));
end
for j=1:1:n-1
u(j)=1-r(j);
end
for j=1:1:n-1
f(j)=(y(j+1)-y(j))/h(j);
end
for j=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);
for j=1:1:n
a(j,j)=2;
end
r(1)=0;
u(n)=0;
for j=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)函数的系数矩阵
for j=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=[0 1 4 9 16 25 36 49 64];
y0=[0 1 2 3 4 5 6 7 8];
x=0:64;y=lagr1(x0,y0,x);
plot(x0,y0,'o')
hold on
plot(x,y,'r');
hold on;
pp=csape(x0,y0,'variational')
fnplt(pp,'g');
hold on;
plot(x0,y0,':b');hold on
%axis([0 2 0 1]); %看[0 1]区间的图形时加上这条指令
gtext('三次样条插值')
gtext('原图像')
gtext('拉格朗日插值')
%下面是求拉格朗日插值的函数
function y=lagr1(x0,y0,x)
n=length(x0); m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=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=[0 0.9 1.9 3.0 3.9 5.0];
y=[0 10 30 50 80 110];
r=zeros(2,2);
for i=1:1:6;
r(1,1)=r(1,1)+1;
end
for i=1:1:6
r(1,2)=r(1,2)+x(i);
end
for i=1:1:6
r(2,1)=r(2,1)+x(i);
end
for i=1:1:6
r(2,2)=r(2,2)+x(i)^2;
end
a=zeros(2,1);
for i=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=[0 0.9 1.9 3.0 3.9 5.0];
y=[0 10 30 50 80 110];
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 给定数据的散点图
通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。具有这种特性的曲线很多，我们选取如下三种函数来拟合：
① 多项式 j (x) = a0 + a1x + … + 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文件：
function F=mf(a,t)
F=a(1)*exp(a(2).*t.^(-1));
在命令窗口中敲入如下代码：
t=[5:5:55]
y=[ 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64];
a0=[0.5,0.5];
a=lsqcurvefit('pf',a0,t,y)
回车后可得：
a(1)= 5.5031 ,a(2)= -8.7872
下面的计算问题用matlab编程实现：
x=[5 10 15 20 25 30 35 40 45 50 55 ];
y0=[1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64];
y=log(y0)
y(1)=0
r=zeros(2,2);
for i=1:1:12;
r(1,1)=r(1,1)+1;
end
for i=2:1:12
r(1,2)=r(1,2)+1/x(i);
end
for i=2:1:12
r(2,1)=r(2,1)+1/x(i);
end
for i=2:1:12
r(2,2)=r(2,2)+1/x(i)^2;
end
a=zeros(2,1);
for i=2:1:12
a(1,1)=a(1,1)+y(i);
a(2,1)=a(2,1)+y(i)*(1/x(i)) ;
end
k=r
t=a
d=inv(r);m=d*a
由matlab软件计算得：a= 3.9864,b= -4.8922。所以y(t)= 3.9864e-4.8922/t。
有下面的语句给出拟合后的图形：
x=[0 5 10 15 20 25 30 35 40 45 50 55 ];
y0=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64];
fplot('5.5031*exp((-8.7872)*(x).^(-1))',[0,55]);
hold on
plot(x,y0,'o',x,y0,'r')
hold on
gtext('倒指数拟合曲线')
gtext('原图像')

图18-2