实验报告数值分析.docx

实验报告数值分析.docx

《实验报告数值分析.docx》由会员分享，可在线阅读，更多相关《实验报告数值分析.docx（50页珍藏版）》请在冰豆网上搜索。

实验报告数值分析

《数值分析》实验报告

姓名：

学号：

专业：

指导教师：

刘建生教授

日期：

2015年12月25日

实验一Lagrange/newton插值

一：

对于给定的一元函数yf（x）的n+1个节点值yjf（Xj）,j0,1,卅，n。

试用Lagrange

公式求其插值多项式或分段二次Lagrange插值多项式。

数据如下：

求

五Xj

欠

yj

计算f（0.596）,f（0.99）的值（提示：

结果为f（0.596）0.625732

f（0.99）1.05423）

Xj

1

2

3

4

5

6

7

yj

试构造Lagrange多项式Lg（x），计算的f（1.8）,f（6.15）值。

（提示：

结果为

f（1.8）0.164762,f（6.15）0.001266）

二:

实验程序及注释

MATLAB程序：

functionf=lagrange（x0,y0,x）

n=length（xO）;

m=length（yO）;

formatlong

s=;

fork=1:

n

P=;

forj=1:

n

ifj~=k

p=p*（x-x0（j））/（x0（k）-x0（j））;

end

end

s=s+y0（k）*p;

End

f=s;

end

结果运行:

>>[O.40-bts0.eb0*SO0.1.Ob]:

yO=L0-410750.57filG0.S067B0.901.00t.253S2J:

x=0»596:

】^Er^nepCvt）,y

ana■

0.6257323?

7352S6a

»xQ=L0,4Q,550,65Q,SO0-951,Q5J;

jrO=TO・410750・578150・696750.90L-00I.253321:

）r=C.99:

lagrari£etarO,yO,vt）

an5—

I、06422&77CS127TE

»xO=ri2345S；]：

»y（J=[0.363C_1350_OSO0.OLE0.007fl.0020.001]:

»x=l.S:

»Lagrangeyfljx）

ans=

0.1S47«1894400000

>>lacrarseCxO^yfl,6^15）

ana=

0.001265S2550039L

结果与提示值完全吻合，说明Lagrange插值多项式的精度是很

高的；

f（x）（xx1）（xx2）（xX3）（XX4）（Xx5）

（x0x1）（x0x2）（x0x3）（x0x4）（x0x5）/

（xx0）（xx1）（xx2）（xx3）（xx4）

（x5x0）（x5x1）（x5x2）（x5x3）（x5x4）

同时，若采用三点插值和两点插值的方法，用三点插值的精度更

高。

若同时采用两点插值，选取的节点距离x越近，精度越高。

三：

采用newton插值进行计算

算法程序如下：

formatlong;

x0=[]；

yo=[];

x=；

n=max（size（xO））;

y=yo（i）;

%disp（y）;

s=1;/

dx=yO;

fori=1:

n-1

dxO=dx;

forj=1:

n-i

dx（j）=（dx0a+1）-dx0（j））/（x0（i+j）-x0（j））;

end

df=dx

（1）;s=s*（x-xO（i））;

y=y+s*df;

%%disp（y）;

end

disp（y）

%计算

运行结果：

df=dx；s=s*（x-xOCi））;

y=yH-£+df:

S%disp

end

dispCy）

.025*********

%计算

一卜-Newtonmterpoistion0AccLratevalue

—*—Piecewiseinterpotatran

oooO

koan書A启j.

绘制出曲线图:

与结果相吻合。

所以newton法和Lagrange法的思想是一样的。

Lagrange适合理论分析，但Lagrange法不女口newton法灵活。

Lagrange女口果节点个数改变，算法需要重新编写，而Newton法克服这一缺点，所以应用更为灵活＜

实验二函数逼近与曲线拟合

一、问题提出

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。

X

t（分钟）0510152025303540455055

/4、

y（10）门

y0

/\

/要求：

1、用最小二乘法进行曲线拟合；\

23

2、近似解析表达式为f（x）=a11+a2t+a3t；

3、计算出拟合函数f（x），并列出出f（x）与y（x）的误差；

4、另外选取一个近似表达式，尝试拟合效果的比较；

5、绘制出曲线拟合图。

二、问题分析

三＞

从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大

量存在，通常利用数据的最小二乘法求得拟合曲线。

三、实验程序及注释

\三次拟合程序（最小二乘法）：

t=[0510152025303540455055]%输入时间t的数据

y=[0]%输入含碳量数据

[p,s]=polyfit（t,y,3）%调用MATLAB最小二乘法的程序进行三次拟合并给出误差分析

formatIong%14位精度小数

plot（t,y,'*r'）%绘制被拟合数据点的离散图

holdon

plot（t,y1,'b'）%绘制三次拟合函数图（其中y1是拟合之后的数据）

xlabel（'时间t（分钟）'）%注释x轴

ylabel（'含碳量/10A-4'）%注释y轴

title（'三次拟合图'）%注释图名

grid%坐标系网格化

四次拟合程序（最小二乘法）：

[p,s]=polyfit（t,y,4）%调用MATLAB最小二乘法的程序进行四次拟合并给出误差分析

formatIong%14位精度小数

plot（t,y,'*r'）%绘制被拟合数据点的离散图

holdon

plot（t,y2,'b'）%绘制三次拟合函数图（其中y2是拟合之后的数据）

xlabel（'时间t（分钟）'）%注释x轴

ylabel（'含碳量/10A-4'）%注释y轴

title（'四次拟合图'）%注释图名

grid%坐标系网格化

四、实验数据结果及分析

三次拟合可以得到其拟合多项式为：

32

y1=拟合函数与被拟合函数图之间的对比如下：

（1）红色星号为原始数据；

（2）带圈的曲线为最小二乘后而成的结果曲线。

由此可见拟合函数与原函数离散数据点拟合成程度相当好，通过[p,s]=polyfit（t,y,n）对拟

合误差进行分析，如图：

图2-2

由此可知，三次拟合精度较好。

为了提高结果的可信度，我们另外选取一个近似表达式，尝试拟合效果的比较。

于是,进行四次拟合：

其中，拟合得到的多项式为：

4

y2=++

拟合如图2-3

7Fi需u"l匚叵|区|

554

4

5352515了21,0.

O

四次拟合图

图2-3

同样对四次拟合进行误差分析可得：

图2-4

由此可见，四次拟合误差<（三次拟合误差），精度更高。

五、实验结论

在用高阶多项式对某一函数进行曲线拟合时，?

并不是拟合出来的多项式与被拟合函数

在整个区间上都能符合，polyfit（）只能保证在输入数据x所能达到的区间上及其附近•求得的

多项式可以最大限度在逼近原函数。

利用最小二乘法对本问题进行的曲线拟合精度较高，而

且，在一般情况下，拟合的多项式次数越多，精度越高。

实验三数值积分与数值微分

、问题提出

（1）

1.5343916）

（2）

（3）

1.

sinx,

dx（f（0）oX

1,10.9460831）

1x

e

I=2dx

o4x

计算下列积分值:

1/4

4-sin2xdx（l

0

1

（4）

.ln（1x）」

I=厂dx

o1x

要求：

1、编制数值积分算法的程序；

2、分别用两种算法计算同一个积分，并比较其结果；

3、分别取不同步长h（ba）/n，试比较计算结果（如n=10,20等）；

4、给定精度要求£,试用变步长算法，确定最佳步长。

/二、问题分析

由上可知这四个积分找不到用初等函数表示的原函数，直接计算起来很困难，因此我们

考虑利用函数在若干点得函数值，近似地计算该函数在一个区间上得定积分。

这里采用的方法有三种：

复合梯形公式，复合Simpson公式，Romberg算法。

三、实验程序及注释

1、复合梯形公式MATLAB程序:

functionl=T_quad（x,y）

%复化梯形求积公式，其中，

%x为向量，被积函数自变量的等距节点；

%y为向量，被积函数节点处的函数值；

n=length（x）;m=length（y）;

ifn〜=m

error（'thelengthofXandYmustbeequal!

'）;return;

end

h=（x（n）-x

（1））/（n-1）;

a=[12*ones（1,n-2）1];

I=h/2*sum（a.*y）;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2、复合Simpson公式MATLAB程序：

functionl=S_quad（x,y）

%x为向量，被积函数自变量的等距节点；

%y为向量，被积函数节点处的函数值；

n=length（x）;m=length（y）;

ifn〜=m

error（'thelengthofXandYmustbeequal!

'）;

return;

end

ifrem（n-1,2）~=0%如果n-1不能被2整除，则调用复化梯形公式

I=T_quad（x,y）;

return;

end

N=（n-1）/2;h=（x（n）-x

（1））/N;a=zeros（1,n）;

fork=1:

N

a（2*k-1）=a（2*k-1）+1;

a（2*k）=a（2*k）+4;

a（2*k+1）=a（2*k+1）+1;

end

l=h/6*sum（a.*y）;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3、Romberg算法MATLAB程序：

functionl=R_quad_iter（fun,a,b,ep）

%Romberg求积公式，其中，

%fun为被积函数；

%a,b为积分区间端点，要求a

%ep精度系数，缺省值为1e-5.

ifnargin<4

ep=1e-5;\/

end

m=1;h=b-a;

l=h/2*（feval（fun,a）+feval（fun