数值分析上机实验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);

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 多项式由下式定义

x(t)??Ct(1?t)kkmk?0mm?kkkxk， y(t)??Cmt(1?t)m?kyk

k?0m

本文来源：https://www.bwwdw.com/article/t54f.html