数值计算(分析)实验报告2

南昌航空大学数学与信息科学学院实验报告

课程名称：《数值计算方法》

实验名称：曲线拟合

实验类型：验证性■综合性□设计性□

实验室名称：数学实验室

班级学号：09072113

学生姓名：邢宪平

任课教师（教师签名）：

成绩：

实验日期：2012/3/25

一、实验目的

实验目的：实验目的：了解函数逼近与曲线拟合的基本原理，并且运用MA TLAB 软件进行实践操作。

二、实验原理、程序框图、程序代码等 实验题目：

题目1：试分别用抛物线2y a bx cx =++和指数曲线bx y ae =拟合下列数据

i x 1 1.5 2 2.5 3 3.5 4 4.5

i y 33.4 79.50 122.65 159.05 189.15 214.15 238.65 252.50 i x 5

5.5

6

6.5

7

7.5

8

i y

267.55 280.50 296.65 301.40 310.40 318.15 325.15

并比较两个拟合函数的优劣。

题目2：已知实验数据如下：

i x 1.0 2.5 3.5 4.0

i y

3.8 1.50 26.0 33.0

试用形如2y a bx =+的抛物线进行最小二乘拟合。

实验原理：

1、逼近方式

假设()[,]f x C a b ∈，2{1,,,...,}n

n H span x x x =，()n n P x H ∈，称

(,)|||||

()()|max n n n a x b f P F P f x P x ≤≤=-=

- 为()f x 与()|n P x 在[,]a b 上的偏差。若存在

*

()n n P x H ∈，使得*

*

(,)|||||()()|max

inf n n

n n n P H a x b

f P f P f x P x ∞∈≤≤=-=

- 则称*

()n P x 是()f x 在[,]a b 上的最

佳一致逼近多项式。

假设()[,]f x C a b ∈及[,]C a b 的一个子集01{(),(),,...()}n span x x x ϕ=ϕϕϕ，若存在*

()S x ϕ∈，使

*

22

2

2

2

()

()||()()||||()()||()[()()]m in m in b

S x S x a

f x S x f x S x x f x S x dx ϕ

ϕρ∈∈-=

-=

-⎰

则称*

()S x 是()f x 在子集[,]

C a b ϕ⊂中的最佳平方逼近数。

2、曲线拟合

上述函数的最佳平方逼近法中，若()f x 是以一组离散点集的形式给出的，即给出了函数()f x 在一些离散点上的值{(,),0,1,...,}i i x y i m =，则该方法就是所说的曲线拟合。取

22

2

(,)||||[()]

m

n n i

i f S f S S x yi ==-=

-∑ 而0011()()()...()n n S x a x a x a x ϕϕϕ=+++，这就是最小二乘

逼近。

由上可知，最小二乘逼近即求多元函数2012

(,,...)||||n n I a a a f S =-的最小值问题。由取最小值的必要条件，可得法方程01()

()...()(,)()()...()(,)

............

......()

()

...

()(,)n n n n n n a y a y n a y 0010000011110ϕ,ϕϕ,ϕϕ,ϕϕϕ,ϕ1ϕ,ϕϕ,ϕϕϕ,ϕϕ1,ϕϕ,ϕϕ⎛⎫⎛⎫⎛⎫

⎪ ⎪ ⎪

⎪ ⎪ ⎪= ⎪ ⎪ ⎪

⎪ ⎪ ⎪

⎝⎭⎝⎭⎝⎭

实际计算时，为了方便，通常取01{(),(),,...()}n span x x x ϕ=ϕϕϕ为简单多项式，如

2

{1,,,...,}n

span x x x ，得到n 次拟合多项式，但当n 较大时，方程组往往是病态的，这时我们可以选

择正交多项式。

程序代码：

题目一

（1）用抛物线y=a+bx+cx^2拟合的程序代码 x=1:0.5:8;

y=[33.4 79.50 122.65 159.05 189.15 214.15 238.65 252.50 267.55 280.50 296.65 301.40 310.40 318.15 325.15]; A=polyfit(x,y,2);

z=polyval(A,x); plot(x,y,'k+',x,z,'r')

（2）用指数曲线y=a*exp{bx}拟合的程序代码

M-文件代码：function f=curvefun1(x,tdata)

f=x(1)*exp{x(2)*tdata} %x(1)=a;x(2)=b;

输入命令： xdata=1:0.5:8

ydata=1e-03*[33.4,79.50,122.65,159.05,189.15,214.15,238,65,252.50,267.55,280.50,296.65,301.40,310.40,318.15,325.15];

t0=[0.2,0.05,0.05];

t=1sqcurvefit('curvefun1',t0,xdata,ydata) f=curvefun1(t,xdata) 题目二程序代码：

function f=curvefun1(x,tdata)

f=x(1)+x(2)*tdata.^2 %其中x （1）=a; x(2)=b;x(3)=k;