数值分析实验

实验三 函数逼近与快速傅里叶变换 P95

专业班级：信计131班 姓名：段雨博 学号：07 一、实验目的

1、熟悉matlab 编程。

2、学习最小二乘法及程序设计算法。 二、实验题目 1、对于给函数()2

1

125f x x

=

+在区间[]1,1-上取()10.20,1,10i x i i =-+=，试求3次曲线拟合，试画出拟合曲线并打印出方程，与第二章计算实习题2的结果进行对比。 x y 图示数据曲线及相应的三种拟合曲线。 3. 给定数据点(),i i x y 如表所示

i x 0 1

i

y

1

用最小二乘法求拟合数据的二次多项式，并求平方误差。 三、实验原理与理论基础

1．最小二乘原理与线性拟合：

在函数的最佳平方逼近中],[)(b a C x f ∈，如果)(x f 只在一组离散点集}...,1,0,{m i x i =上给出，这就是科学实验中经常见到的实验数据｛}...1,0),,{(m i y x i i =｝的曲线拟合，这里

)...1,0)((m i x f y i i ==，要求一个函数)(*x S y =与所给数据}...1,0),,{(m i y x i i =拟合，

若记误差T m ),...(10

δδδδ=，设)(),...(),(10x x x n ϕϕϕ是C[a,b]上线性无关函数族，在

)}(),...(),({10x x x span n ϕϕϕϕ=中找一函数)(*x S 使误差平方和

2

20

222])([min ])(*[∑∑∑===-=-==m

i i i m i i m i i y x S y x S δδ

，这

里 )(...)()()(1000x a x a x a x S n n ϕϕϕ+++= （m n <）。这就是一般的最小二乘逼近，用几何语言说，就称为曲线拟合的最小二乘法。

数据拟合是根据测定的数据间的相互关系，确定曲线),...,,;(10n a a a x s y =的类型，然后再根据在给定点上误差的平方和达到最小的原则，即求解无约束问题：

∑=-=m

i i n i n y a a a x s a a a F 1

21010,)),...,,;((),...,,(m in

确定出最优参数：),...,1,0(*n k a k

=，从而得到拟合曲线).(*

x s y =

2、多项式拟合。

3、定义1：

设有数据},...,,{21m x x x X =和权系数),,...,2,1(m i i =ω称：

∑==

m

i i i i

x g x f g f 1

)()(),(ω

为函数为权的内积。上以在和ω},...,,{21m x x x X g f =

4、用正交函数最佳平方逼近：

为避免出现正规方程组的系数矩阵是病态矩阵的情况，在选择多项式时需要考虑正交的多项式n n b a C span φφφφφφ,...,,],,[},...,,{.1010⊂=Φ设是正交基，即：

⎪⎩

⎪⎨⎧=≠=k i k

i k k i ,,0),(22φφφ 于是正规方程组可以简化为：,,...,1,0),,(2

2n k f a k k k ==φφ （1）

解方程得到：.,...,1,0,/),(2

2*

n k f a a k k k k ===φφ这里避免了求解病态方程组，提高了计

算系数的精确度。对任意的],,[)(b a C x f ⊂其最佳平方逼近函数为：

)(/),()()(22

*

*

x f x a x s k k

n

k n

k k k k

φφφφ∑∑===

由式（1）以及∑=-=-+-=-=n

k k k f a f

s f s f s f

s

f 0

*

22

*

*

*2

2

*22

),(),(),(φδ，导出平方误差为：

∑=-=n

k k k

a f

2

2*2222

,)(φδ

其平方根称为均方误差。

5、由题意决定,....),,1(2

x x span ,即决定拟合多项式，分别计算1

(,)(,)n

i j i

j

m k k k k ==

∑，

1

(,)(,)n i i j m k y k y -

==∑，用(,)i j k k 组成方阵A,用(,)i k y -

组成矩阵B ，利用A/B 求出该多项式

的系数，再利用

1

(())^2n

i

i

i f x y =-∑求出平方误差。

四、实验内容

解：1、 >> i = 0:10; >> x = -1+*i;

>> y = 1./(1+25*x.^2); >> p=polyfit(x,y,3); >> s=vpa(poly2sym(p)) s =

- *x^3 - *x^2 + *x + >> f=polyval(p,x); >> plot(x,f,x,y,'o ')

2、

>> x=[0 1];

>> y=[1 ];

>> p1=polyfit(x,y,3)%三次多项式拟合

p1 =

>> p2=polyfit(x,y,4)%四次多项式拟合

p2 =

>> y1=polyval(p1,x);

>> y2=polyval(p2,x);%多项式求值

>> plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')

>> p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。p3 =

>> y3=polyval(p3,x);

>> plot(x,y,'c--',x,y1,'r:',x,y2,'y-.',x,y3,'k--')%画出四种拟合曲线

3、M文件：

function []=zuixiaoercinihe2(x,y)

n=length(x);

k00=0;

for i=1:n

k00=k00+1;

end

k01=0;

for i=1:n

k01=k01+x(i);

end

k02=0;

for i=1:n

k02=k02+x(i)*x(i);

end

k11=0;

for i=1:n

k11=k11+x(i)*x(i);

end

k12=0;

for i=1:n

k12=k12+x(i)*x(i)*x(i);

end

k22=0;

for i=1:n

k22=k22+x(i)*x(i)*x(i)*x(i);