数值分析上机实验6

1、数值分析上机实验6例1世界人口数据拟合问题：据统计，六十年代世界人口数据如下(单位：亿)年196019611962196319641965196619671968人口29.7230.6131.5132.1332.3432.8533.5634.2034.83根据表中数据，预测公元2000年时的世界人口。问题分析与数学模型设人口总数为 N(t)，根据人口理论的马尔萨斯模型， 采用指数函数N(t) = e a + b t对数据进行拟合。为了计算方便，将上式两边同取对数，得 ，令 y = ln N 或 N = e y变换后的拟合函数为y(t) = a + b t由人口数据取对数（y = ln N）计算

2、，得下表t196019611962196319641965196619671968y3.39183.42133.45033.46983.47633.49203.51333.53223.5505根据表中数据及等式a + b t k = y k ( k = 1，2，9)可列出关于两个未知数a 、b的9个方程的超定方程组（方程数多于未知数个数的方程组）a + tj b = yj (j= 1，2，9)可用最小二乘法求解。算法与数学模型求解算法如下：第一步：输入人口数据，并计算所有人口数据的对数值；第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；第三步：求解超定方程组并输出

3、结果：a，b；第四步：利用数据结果构造指数函数计算2000年人口近似值N(2000)，结束。MATLAB程序t=1960:1968;t0=2000;N=29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.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.0186N(2000)

4、= 63.2336所以取五位有效数，可得人口数据的指数拟合函数经计算得2000年人口预测值为：63.2336 (亿)。例2温度数据的三角函数拟合问题洛杉矶郊区在11月8日的温度记录如下时间1234567891011午夜温度666665646363626160605958时间1314151617181920212223正午温度585858585757575860646768在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier级数的部分和（有限项）做拟合函数。即，求最小二乘曲线：其中，w = 2 / 24。例如，当n=2时拟合函数为j(x)= a0 + a1 cos(wx) + b1

5、 sin(wx) + a2 cos(2wx) + b2 sin(2wx)对不同的n，确定拟合函数中的各系数。绘出最小二乘曲线与离散数据点，并计算出拟合函数的残差2-范数。算法分析：以n=2时的拟合函数为对象作算法分析。24小时的温度记录可列为数表如下xx1 x2 x3 x24yy1 y1 y1 y24将24个数据点代入拟合函数得超定方程组可以证明方程组的系数矩阵列向量组是正交向量组，于是由最二乘法所推出的正规方程组系数矩阵是对角矩阵。所以原方程组的最小二乘解为， MATLAB程序（运行程序时需输入参数n ）：n=input('input n=: ');w=2*pi/24;x=1

6、: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;for k=1:n ck=cos(k*w*x);sk=sin(k*w*x); a(k)=(ck'*y)/(ck'*ck); b(k)=(sk'*y)/(sk'*sk);endyy=a0;for k=1:n yy=yy+a(k)*cos(k*w*x)+b(k)*sin(k*w*x);endplot(x,y,'x',x,yy)r=norm(yy-y)n=1，r=

7、 7.32853切比雪夫多项式的前两项为：T0(x) = 1，T1(x) = x，对于n2，有递推公式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');hold onfor k=2:4 T=2*x.*T1-T0; plot(x,T) T0=T1;T1=T;endaxis off41912年，伯恩斯坦给出了关于多项式一致逼

8、近连续函数的构造性证明，提出了著名的伯恩斯坦多项式，设 f(x)在区间 0，1上连续，他的多项式为试利用组合数的递推公式 ，设计一个计算n次伯恩斯坦多项式函数值的算法。并对函数 f(x) = sin x 给以验证。MATLAB程序如下n=input('input n=');x=0:n/n; f=sin(x*pi);for i=1:n+1 y=f;t=x(i); for k=n:-1:1 for j=1:k y(j)=t*y(j)+(1-t)*y(j+1); end end p(i)=y(1);endmax(abs(f-p)plot(x,f,'b',x,p,

9、9;o',x,p,'r')运行四次程序，分别输入 n=10，n=20，n=30，n=40得下面图形 5Bezier 曲线是法国雷诺汽车公司的工程师Bezier 于1971年提出了一种新的参数曲线表示法。这种方法可以交互式地确定一组控制多边形顶点以获得所需要的曲线形式。设曲线参数方程x = x(t)， y = y(t)如果给定控制多边形顶点P0，P1，Pm 的坐标(x0，y0)，(x1，y1)，(xm，ym)则相应的Bezier 多项式由下式定义， 用矢量函数的形式来表示平面曲线则有数学表达式三次Bezier 曲线数学表达式为 试给定四个平面点绘制飞机机翼剖面轮廓曲线fu

10、nction z=bez3(p)t=0:.05:1't1=1-t;z=t1.3 3*t.*t1.2 3*t1.*t.2 t.3*p;plot(p(:,1),p(:,2),z(:,1),z(:,2)首先录入这一函数文件，文件名为bez3.m。为了调用这一函数，必须先定义多边折线的结点数组，在MATLAB环境下键入：p=0 0;0 0.4;0.5 1;2 0bez3(p)计算机将绘制出对应的曲线。一、数值试验1. SARS的传播及预防问题非典的爆发和蔓延给我国的经济发展和人民生活带来了很大的影响，下表给出了北京市当年4月份到6月份的疫情数据,通过拟合确诊的累积病人曲线,若延后5天采取严格的

11、预防措施,对疫情的传播所生成的影响做出估计.日期已确诊病例累积现有疑似病例死亡累积治愈出院累积4月20日29740218334月30日1584140875905月1日16401415821005月10日198813971161755月20日218912251503955月30日230970617610066月1日231973918111246月10日239435118417476月20日243931912189(1) 根据已有数据，用MATLAB绘出确诊累积病人离散图(2) 取拟合曲线的拟合函数为如下非线性函数试确定拟合函数中的参数：a，b，并推测五天后累积病人数量。2蠓虫分类问题生物学家试图

12、对两类蠓虫(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.08Apf 数据 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）依据确

13、立的方法,对题目提供的三个样本：(1.24,1.80)，(1.28,1.84)，(1.40,2.04)加以识别；（3）设Af是宝贵的传粉益虫，Apf是某种疾病的载体，是否应该修改分类方法。问题分析：首先画出15对数据的散点图，其中，Af 用 标记，Apf 用 × 标记。观察图形，可以发现，Af 的点集中在图中右下角，而 Apf 的点集中在图中左上角。应该存在一条直线 L 位于两类点之间， 作为 Af 和 Apf 分界线，这条直线 L 的确定应依据问题所给的数据，即学习样本。设这条直线的方程为w1 x + w2 y + w0 = 0对于平面上任意一点 P(x，y)，如果该点在直线上，将

14、其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点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)= w1 x +

15、 w2 y + w0其中，w1，w2，w0为三个待定系数，由所给数据形成约束条件，即线性方程组： 这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。实验要求：（1）根据上面分析写出对应的正规方程组并求解；（2）确定分类边界直线的方程。由所给数据用判别函数判别三个新蠓虫的类属；（3）如果Apf是某种疾病的载体，应该修改数学模型确定分类方法。3海浪高度数据的拟合问题海洋水文观测站记录了某海域每天24小时海浪潮高度数据（相对于海堤上的零标尺记号）。工作人员每小时做一次记录，下面是从某年12月1日早上零点开始到12月2日晚上23点结束两天内完整的数据记录。能否根据这些数据预测三天以

16、后（即12月5日）的海浪高度数据。182.40001.2000-0.1000-0.5000-2.5000-3.0000-2.7000-1.60009160.20002.10003.40003.60002.90001.60000.2000-1.20001724-2.4000-3.0000-3.1000-2.3000-0.70001.30002.90003.600025323.10002.00000.60000.6000-2.2000-3.6000-3.2000-2.50003340-0.9000-1.10002.90003.90003.60002.50001.0000-1.50004148-2.

17、4000-3.0000-3.4000-3.0000-1.70000.20002.20003.5000 二、方法应用与编程练习1龙格（Range）反例是插值方法中一个著名的反例。比较多项式插值和多项式拟合的方法的异同，借用MATLAB中数据拟合命令做龙格反例的数值实验。2人口预测问题。已知我国人口数据如下（单位：亿）年t199119921993199419951996数量N11.5811.7211.8511.9812.1112.24根据马尔萨斯人口模型，设N(t) = exp( a + b t )，利用数据拟合方法确定指数函数，并预测2000年到2008年九年的我国人口数量。3勒让德多项式的前两

18、项为：p0(x) = 1，p1(x) = x，对于n2，有递推公式pn+1(x) = (2n+1) x pn(x) n pn 1(x)/(n+1)对x 1，1 ，绘制勒让德多项式p0(x)，p1(x)，p2(x)，p3(x)，p4(x)图形4假设慧星1968Tenta在太阳系内移动，在某个极坐标系下对慧星位置记录下观察值如下r2.702.001.611.201.02j48°67°83°108°126°由Kepler第一定律，慧星应在一个椭圆或双曲线的平面轨道上运行，忽略干扰有轨道方程其中，p 为参数，e为偏心率，由给定的观察值用最小二乘方法拟合轨道曲线，求出p和e，并给出平方误差。 令k=1/r，处理数据。解超定方程组确定p和er=2.70 2.00 1.61 1.20 1.02;fai=48 67 83 108 126*pi/180;n=length(r);b=ones(n,1);G=1./r' cos(fai');X=Gb;p=X(1);e=X(2);fa=pi/4:.1:3.5*pi/2;ro=p./(1-e*cos(fa);polar(fai,r,'ro'),hold onpolar(fa,ro)5斐波拉奇数列中整