上机实习六.docx
《上机实习六.docx》由会员分享，可在线阅读，更多相关《上机实习六.docx（15页珍藏版）》请在冰豆网上搜索。
上机实习六
数值分析上机实验6
例1．世界人口数据拟合问题：
据统计，六十年代世界人口数据如下（单位：
亿）
年
1960
1961
1962
1963
1964
1965
1966
1967
1968
人口
29.72
30.61
31.51
32.13
32.34
32.85
33.56
34.20
34.83
根据表中数据，预测公元2000年时的世界人口。
问题分析与数学模型
设人口总数为N（t），根据人口理论的马尔萨斯模型，采用指数函数
N（t）=ea+bt
对数据进行拟合。
为了计算方便，将上式两边同取对数，得
，令
y=lnN或N=ey
变换后的拟合函数为
y（t）=a+bt
由人口数据取对数（y=lnN）计算，得下表
t
1960
1961
1962
1963
1964
1965
1966
1967
1968
y
3.3918
3.4213
3.4503
3.4698
3.4763
3.4920
3.5133
3.5322
3.5505
根据表中数据及等式a+btk=yk（k=1，2，……，9）可列出关于两个未知数a、b的9个方程的超定方程组（方程数多于未知数个数的方程组）
a+tjb=yj（j=1，2，…，9）
可用最小二乘法求解。
算法与数学模型求解
算法如下：
第一步：
输入人口数据，并计算所有人口数据的对数值；
第二步：
建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；
第三步：
求解超定方程组并输出结果：
a，b；
第四步：
利用数据结果构造指数函数计算2000年人口近似值N（2000），结束。
MATLAB程序
t=1960:
1968;t0=2000;
N=[29.7230.6131.5132.1332.3432.8533.5634.2034.83];
y=log（N）;
A=[ones（9,1）,t'];
d=A\y';a=d
（1）,b=d
（2）
N0=exp（a+b*t0）
x=1960:
2001;yy=exp（a+b*x）;
plot（x,yy,t,N,'o',2000,N0,'o'）
计算结果为
－33.0383，
0.0186
N（2000）=63.2336
所以取五位有效数，可得人口数据的指数拟合函数
经计算得2000年人口预测值为：
63.2336（亿）。
例2．温度数据的三角函数拟合问题
洛杉矶郊区在11月8日的温度记录如下
时间
1
2
3
4
5
6
7
8
9
10
11
午夜
温度
66
66
65
64
63
63
62
61
60
60
59
58
时间
13
14
15
16
17
18
19
20
21
22
23
正午
温度
58
58
58
58
57
57
57
58
60
64
67
68
在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier级数的部分和（有限项）做拟合函数。
即，求最小二乘曲线：
其中，=2π/24。
例如，当n=2时拟合函数为
（x）=a0+a1cos（x）+b1sin（x）+a2cos（2x）+b2sin（2x）
对不同的n，确定拟合函数中的各系数。
绘出最小二乘曲线与离散数据点，并计算出拟合函数的残差2-范数。
算法分析：
以n=2时的拟合函数为对象作算法分析。
24小时的温度记录可列为数表如下
x
x1x2x3…………x24
y
y1y1y1…………y24
将24个数据点代入拟合函数得超定方程组
可以证明方程组的系数矩阵列向量组是正交向量组，于是由最二乘法所推出的正规方程组系数矩阵是对角矩阵。
所以原方程组的最小二乘解为
，
，
MATLAB程序（运行程序时需输入参数n）：
n=input（'inputn=:
'）;
w=2*pi/24;x=[1:
24]';
y=[66;66;65;64;63;63;62;61;60;60;59;58;58;58;58;58;57;57;57;58;60;64;67;68];
a0=sum（y）/24;
fork=1:
n
ck=cos（k*w*x）;sk=sin（k*w*x）;
a（k）=（ck'*y）/（ck'*ck）;b（k）=（sk'*y）/（sk'*sk）;
end
yy=a0;
fork=1:
n
yy=yy+a（k）*cos（k*w*x）+b（k）*sin（k*w*x）;
end
plot（x,y,'x',x,yy）
r=norm（yy-y）
n=1，r=7.3285
3．切比雪夫多项式的前两项为：
T0（x）=1，T1（x）=x，对于n≥2，有递推公式
Tn+1（x）=2xTn（x）–Tn–1（x）
当x∈[–1，1]时，利用递推公式，计算并绘出T0（x），T1（x），T2（x），T3（x），T4（x）的函数图形
MATLAB程序如下：
x=-1:
.05:
1;
T0=ones（size（x））;
T1=x;
plot（x,T0,'b',x,T1,'b'）;
holdon
fork=2:
4
T=2*x.*T1-T0;
plot（x,T）
T0=T1;T1=T;
end
axisoff
4．1912年，伯恩斯坦给出了关于多项式一致逼近连续函数的构造性证明，提出了著名的伯恩斯坦多项式，设f（x）在区间[0，1]上连续，他的多项式为
试利用组合数的递推公式
，设计一个计算n次伯恩斯坦多项式函数值的算法。
并对函数f（x）=sinx给以验证。
MATLAB程序如下
n=input（'inputn='）;
x=[0:
n]/n;f=sin（x*pi）;
fori=1:
n+1
y=f;t=x（i）;
fork=n:
-1:
1
forj=1:
k
y（j）=t*y（j）+（1-t）*y（j+1）;
end
end
p（i）=y
（1）;
end
max（abs（f-p））
plot（x,f,'b',x,p,'o',x,p,'r'）
运行四次程序，分别输入n=10，n=20，n=30，n=40得下面图形
5．Bezier曲线是法国雷诺汽车公司的工程师Bezier于1971年提出了一种新的参数曲线表示法。
这种方法可以交互式地确定一组控制多边形顶点以获得所需要的曲线形式。
设曲线参数方程
x=x（t），y=y（t）
如果给定控制多边形顶点P0，P1，…，Pm的坐标
（x0，y0），（x1，y1），……，（xm，ym）
则相应的Bezier多项式由下式定义
，
用矢量函数的形式来表示平面曲线则有数学表达式
三次Bezier曲线数学表达式为
试给定四个平面点绘制飞机机翼剖面轮廓曲线
functionz=bez3（p）
t=[0:
.05:
1]';t1=1-t;
z=[t1.^33*t.*t1.^23*t1.*t.^2t.^3]*p;
plot（p（:
1）,p（:
2）,z（:
1）,z（:
2））
首先录入这一函数文件，文件名为bez3.m。
为了调用这一函数，必须先定义多边折线的结点数组，在MATLAB环境下键入：
p=[00;00.4;0.51;20]
bez3（p）
计算机将绘制出对应的曲线。
一、数值试验
1.SARS的传播及预防问题
非典的爆发和蔓延给我国的经济发展和人民生活带来了很大的影响，下表给出了北京市当年4月份到6月份的疫情数据,通过拟合确诊的累积病人曲线,若延后5天采取严格的预防措施,对疫情的传播所生成的影响做出估计.
日期
已确诊病例累积
现有疑似病例
死亡累积
治愈出院累积
4月20日
297
402
18
33
4月30日
1584
1408
75
90
5月1日
1640
1415
82
100
5月10日
1988
1397
116
175
5月20日
2189
1225
150
395
5月30日
2309
706
176
1006
6月1日
2319
739
181
1124
6月10日
2394
351
184
1747
6月20日
2439
3
191
2189
（1）根据已有数据，用MATLAB绘出确诊累积病人离散图
（2）取拟合曲线的拟合函数为如下非线性函数
试确定拟合函数中的参数：
a，b，并推测五天后累积病人数量。
2．蠓虫分类问题
生物学家试图对两类蠓虫（Af与Apf）进行鉴别,依据的资料是蠓虫的触角和翅膀的长度,已经测得9只Af和6只Apf的数据（触角长度用x表示,翅膀长度用y表示）
Af数据
x
1.24
1.36
1.38
1.38
1.38
1.40
1.48
1.54
1.56
y
1.27
1.74
1.64
1.82
1.90
1.70
1.82
1.82
2.08
Apf数据
x
1.14
1.18
1.20
1.26
1.28
1.30
y
1.78
1.96
1.86
2.00
2.00
1.96
现需要解决三个问题：
（1）如何凭借原始资料（15对数据,被称之为学习样本）制定一种方法区分两类蠓虫；
（2）依据确立的方法,对题目提供的三个样本：
（1.24,1.80），（1.28,1.84），（1.40,2.04）加以识别；
（3）设Af是宝贵的传粉益虫，Apf是某种疾病的载体，是否应该修改分类方法。
问题分析：
首先画出15对数据的散点图，其中，Af用
标记，Apf用×标记。
观察图形，可以发现，Af的点集中在图中右下角，而Apf的点集中在图中左上角。
应该存在一条直线L位于两类点之间，作为Af和Apf分界线，这条直线L的确定应依据问题所给的数据，即学习样本。
设这条直线的方程为
w1x+w2y+w0=0
对于平面上任意一点P（x，y），如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点P（x，y）不在直线上，将其坐标代入直线方程，则方程左端不为零。
由于Af和Apf的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果。
为了建立判别准则，利用直线方程左端表达式定义判别函数g（x，y），使得当平面上某点P（x，y）属于Af类时，g（x，y）>0，否则g（x，y）<0。
引入学习机制，在学习过程中将两种不同的状态，以“1”和“－1”表示。
当P（x，y）属于Af类时，g（x，y）=1，否则g（x，y）=－1。
取
g（x，y）=w1x+w2y+w0
其中，w1，w2，w0为三个待定系数，由所给数据形成约束条件，即线性方程组：
这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。
实验要求：
（1）根据上面分析写出对应的正规方程组并求解；
（2）确定分类边界直线的方程。
由所给数据用判别函数判别三个新蠓虫的类属；
（3）如果Apf是某种疾病的载体，应该修改数学模型确定分类方法。
3．海浪高度数据的拟合问题
海洋水文观测站记录了某海域每天24小时海浪潮高度数据（相对于海堤上的零标尺记号）。
工作人员每小时做一次记录，下面是从某年12月1日早上零点开始到12月2日晚上23点结束两天内完整的数据记录。
能否根据这些数据预测三天以后（即12月5日）的海浪高度数据。
1~8
2.4000
1.2000
-0.1000
-0.5000
-2.5000
-3.0000
-2.7000
-1.6000
9~16
0.2000
2.1000
3.4000
3.6000
2.9000
1.6000
0.2000
-1.2000
17~24
-2.4000
-3.0000
-3.1000
-2.3000
-0.7000
1.3000
2.9000
3.6000
25~32
3.1000
2.0000
0.6000
0.6000
-2.2000
-3.6000
-3.2000
-2.5000
33~40
-0.9000
-1.1000
2.9000
3.9000
3.6000
2.5000
1.0000
-1.5000
41~48
-2.4000
-3.0000
-3.4000
-3.0000
-1.7000
0.2000
2.2000
3.5000
二、方法应用与编程练习
1．龙格（Range）反例是插值方法中一个著名的反例。
比较多项式插值和多项式拟合的方法的异同，借用MATLAB中数据拟合命令做龙格反例的数值实验。
2．人口预测问题。
已知我国人口数据如下（单位：
亿）
年t
1991
1992
1993
1994
1995
1996
数量N
11.58
11.72
11.85
11.98
12.11
12.24
根据马尔萨斯人口模型，设N（t）=exp（a+bt），利用数据拟合方法确定指数函数，并预测2000年到2008年九年的我国人口数量。
3．勒让德多项式的前两项为：
p0（x）=1，p1（x）=x，对于n≥2，有递推公式
pn+1（x）=[（2n+1）xpn（x）–npn–1（x）]/（n+1）
对x∈[–1，1]，绘制勒让德多项式p0（x），p1（x），p2（x），p3（x），p4（x）图形
4．假设慧星1968Tenta在太阳系内移动，在某个极坐标系下对慧星位置记录下观察值如下
r
2.70
2.00
1.61
1.20
1.02
48
67
83
108
126
由Kepler第一定律，慧星应在一个椭圆或双曲线的平面轨道上运行，忽略干扰有轨道方程
其中，p为参数，e为偏心率，由给定的观察值用最小二乘方法拟合轨道曲线，求出p和e，并给出平方误差。
令k=1/r，处理数据。
解超定方程组确定p和e
r=[2.702.001.611.201.02];
fai=[486783108126]*pi/180;
n=length（r）;b=ones（n,1）;
G=[1./r'cos（fai'）];
X=G\b;
p=X
（1）;e=X
（2）;
fa=pi/4:
.1:
3.5*pi/2;
ro=p./（1-e*cos（fa））;
polar（fai,r,'ro'）,holdon
polar（fa,ro）
5．斐波拉奇数列中整数的高精度计算问题。
对n=100和n=200，分别计算斐波拉奇数列的第n项和前n项和（要求全部用整数表示）。
6．多边形面积的计算问题。
任意一个n边形都可以分割成为n-2个三角形，分别计算出这些三角形的面积，将所有三角形面积和作为n边形面积是一种解决问题的有效方法。
输入数据的第一行有一个正整数n，表示多边形顶点数目（0≤n≤1000）。
接下来有n行数据，每一行有一对实数xi，yi用来表示多边形顶点的平面坐标（假定这些顶点确实构成多边形）。
对输入数据计算出多边形面积S。
SampleImput
5
00
01
0.50.5
11
10
OutputfortheSampleInput
0.75
升级会员 