数值分析上机实验6

数值分析上机实验6

例1．世界人口数据拟合问题：据统计，六十年代世界人口数据如下(单位：亿)

年 人口 1960 29.72 1961 30.61 1962 31.51 1963 32.13 1964 32.34 1965 32.85 1966 33.56 1967 34.20 1968 34.83 根据表中数据，预测公元2000年时的世界人口。

问题分析与数学模型

设人口总数为 N(t)，根据人口理论的马尔萨斯模型， 采用指数函数

N(t) = e a + b t

对数据进行拟合。为了计算方便，将上式两边同取对数，得 lnN?a?bt，令

y = ln N 或 N = e y

变换后的拟合函数为

y(t) = a + b t

由人口数据取对数（y = ln N）计算，得下表 t y 1960 1961 1962 1963 1964 1965 1966 1967 1968 3.3918 3.4213 3.4503 3.4698 3.4763 3.4920 3.5133 3.5322 3.5505 根据表中数据及等式a + b t k = y k ( k = 1，2，??，9)可列出关于两个未知数a 、b的9个方程的超定方程组（方程数多于未知数个数的方程组）

a + tj b = yj (j= 1，2，?，9)

可用最小二乘法求解。

算法与数学模型求解 算法如下：

第一步：输入人口数据，并计算所有人口数据的对数值；

第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量； 第三步：求解超定方程组并输出结果：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')

计算结果为

a?－33.0383， b? 0.0186

N(2000) = 63.2336

所以取五位有效数，可得人口数据的指数拟合函数

N(t)?e?33.0383?0.0186t

经计算得2000年人口预测值为：63.2336 (亿)。

70 60 50 40 30 20 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下 时间 1 温度 66 时间 13 温度 58 2 66 14 58 3 65 15 58 4 64 16 58 5 63 17 57 6 63 18 57 7 62 19 57 8 61 20 58 9 60 21 60 10 60 22 64 11 59 23 67 午夜 58 正午 68 在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier级数的部分和（有限项）做拟合函数。即，求最小二乘曲线：

?(x)?a0??[akcos(k?x)?bksin(k?x)]

k?1n其中，? = 2π / 24。例如，当n=2时拟合函数为

?(x)= a0 + a1 cos(?x) + b1 sin(?x) + a2 cos(2?x) + b2 sin(2?x)

对不同的n，确定拟合函数中的各系数。绘出最小二乘曲线与离散数据点，并计算出拟合函

数的残差2-范数。 算法分析：

以n=2时的拟合函数为对象作算法分析。24小时的温度记录可列为数表如下

x y x1 x2 x3 ???? x24 y1 y1 y1 ???? y24 将24个数据点代入拟合函数得超定方程组

?1cos?x1?1cos?x2??????????1cos?x24sin?x1sin?x2??sin?x24cos2?x1cos2?x2??cos2?x24sin2?x1??a0??y1?????aysin2?x2???1??2???b1????? ?????????a2????????sin2?x24???b2??y24?可以证明方程组的系数矩阵列向量组是正交向量组，于是由最二乘法所推出的正规方程组

系数矩阵是对角矩阵。所以原方程组的最小二乘解为

124a0??yk

24k?1a1??cos?xkyk/?cos?xk，a2??cos2?xkyk/?cos22?xk

224242424k?124k?124k?124k?124b1??sin?xkyk/?sin2?xk， b2??sin2?xkyk/?sin22?xk

k?1k?1k?1k?1MATLAB程序（运行程序时需输入参数n ）： n=input('input n=: '); 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; 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); end yy=a0; for k=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)

70 65 60 55 0 5 10 15 20 25 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'); hold on for k=2:4

T=2*x.*T1-T0; plot(x,T)

T0=T1;T1=T; end axis off

4．1912年，伯恩斯坦给出了关于多项式一致逼近连续函数的构造性证明，提出了著名的伯恩斯坦多项式，设 f(x)在区间 [0，1]上连续，他的多项式为

kBn(x)??f()Cnk(1?x)n?kxk

nk?0kkk?1试利用组合数的递推公式 Cn?Cn?1?Cn?1，设计一个计算n次伯恩斯坦多项式函数值的

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);