数值分析上机实验6

数值分析上机实验6

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

根据表中数据，预测公元2000年时的世界人口。

问题分析与数学模型

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

N(t) = e a + b t

=+，令

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

y = ln N或N = e y

变换后的拟合函数为

y(t) = a + b t

根据表中数据及等式 k k( 1，2，……，9)可列出关于两个未知数、b的9个方程的超定方程组（方程数多于未知数个数的方程组）

a + t j b = y j(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

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

t e t N 0186.00383.33)(+-=

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

例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下

在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier 级数的部分和（有限项）做拟合函数。即，求最小二乘曲线：

∑=++=n

k k k x k b x k a a x 1

0)]sin()cos([)(ωωϕ

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

ϕ(x )= a 0 + a 1 cos(ωx ) + b 1 sin(ωx ) + a 2 cos(2ωx ) + b 2 sin(2ωx )

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

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

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

将24个数据点代入拟合函数得超定方程组

⎥⎥

⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢

⎢⎢

⎢⎣⎡2421221102424

2424

22

2

2111

12sin 2cos sin cos 12sin 2cos sin cos 12sin 2cos sin cos 1y y y b a b a a x x x x x x x x

x x x x

ωωωωωωωωωωωω 可以证明方程组的系数矩阵列向量组是正交向量组，于是由最二乘法所推出的正规方程组

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

∑==24

1

0241k k y a

∑∑===24

12

24

11cos /cos k k k k k x y x a ωω，∑∑===24

1224

122cos /2cos k k k k k x y x a ωω

∑∑===24

1

224

1

1sin /sin k k k k k x y x b ωω， ∑∑===24

1

224

1

22sin /2sin k k k k k x y x b ωω

MATLAB 程序（运行程序时需输入参数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)

3．切比雪夫多项式的前两项为：T 0(x ) = 1，T 1(x ) = x ，对于n ≥2，有递推公式

T n+1(x ) = 2xT n (x ) – T n – 1(x )

当x ∈[ – 1，1 ] 时，利用递推公式，计算并绘出 T 0(x )，T 1(x )，T 2(x )，T 3(x )，T 4(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]上连续，他的多项式为

B x f k

n C x x n n k n k k k n

()()()=--=∑10

试利用组合数的递推公式 1

11---+=k n k n k n C C C ，设计一个计算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);