实验报告—数值分析

江西理工大学上数值分析期末实验

《数值分析》实验报告

姓 名： 学 号： 专 业：

指导教师： 刘 建 生 教 授 日 期： 2015年12月25日

江西理工大学上数值分析期末实验

实验一 Lagrange/newton插值

一：对于给定的一元函数y f(x)的n+1个节点值yj f(xj),j 0,1, ,n。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。

数据如下：

6f(0.99)的值（提示：结果为f(0.596) 0.625732, 计算f(0.59,f(0.99) 1.05423 ）

试构造Lagrange多项式L6(x)，计算的f(1.8),f(6.15)值。（提示：结果为

f(1.8) 0.164762, f(6.15) 0.001266 ）

二：实验程序及注释

MATLAB程序：function f=lagrange(x0,y0,x )

n=length(x0); m=length(y0); format long s=0.0; for k=1:n p=1.0; for j=1:n if j~=k

p=p*(x-x0(j))/(x0(k)-x0(j)); end end

s=s+y0(k)*p; End f=s; end

江西理工大学上数值分析期末实验

结果运行：

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

f(x)

(x x1)(x x2)(x x3)(x x4)(x x5)

(x0 x1)(x0 x2)(x0 x3)(x0 x4)(x0 x5)

(x x0)(x x1)(x x2)(x x3)(x x4)(x5 x0)(x5 x1)(x5 x2)(x5 x3)(x5 x4)

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

江西理工大学上数值分析期末实验

三：采用newton插值进行计算 算法程序如下：

format long; x0=[0.4

0.55 0.65 0.80 0.95 1.05 ];

0.69675

0.90 1.00 1.25382 ];

y0=[0.41075 0.57815 n=max(size(x0));

y=y0(1); %disp(y); s=1;

dx=y0; for i=1:n-1 dx0=dx; for j=1:n-i

dx(j)=(dx0(j+1)-dx0(j))/(x0(i+j)-x0(j)); end df=dx(1); s=s*(x-x0(i));

y=y+s*df; %计算 %%disp(y); end disp(y)

x=0.596;

运行结果：

江西理工大学上数值分析期末实验

绘制出曲线图：

与结果相吻合。

所以newton法和Lagrange法的思想是一样的。Lagrange适合理论分析，但Lagrange法不如newton法灵活。Lagrange如果节点个数改变，算法需要重新编写，而Newton法克服这一缺点，所以应用更为灵活。

江西理工大学上数值分析期末实验

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

一、问题提出

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。 t(分钟) y（10 4）

要求：

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

2、近似解析表达式为f(x)=a1t+ a2t2+ a3t；

3、计算出拟合函数f(x)，并列出出f(x)与y(x)的误差； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、绘制出曲线拟合图。

二、问题分析 三、

从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。

三、实验程序及注释

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

t=[0 5 10 15 20 25 30 35 40 45 50 55]%输入时间t的数据

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]%输入含碳量数据

[p,s]=polyfit(t,y,3)%调用MATLAB最小二乘法的程序进行三次拟合并给出误差分析 format long%14位精度小数

plot(t,y,'*r')%绘制被拟合数据点的离散图 hold on

plot(t,y1,'b')%绘制三次拟合函数图（其中y1是拟合之后的数据） xlabel('时间t（分钟）') %注释x轴 ylabel('含碳量/10^-4') %注释y轴 title('三次拟合图') %注释图名 grid%坐标系网格化

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

[p,s]=polyfit(t,y,4) %调用MATLAB最小二乘法的程序进行四次拟合并给出误差分析 format long%14位精度小数

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

3

0 0

5 1.27

10 2.16

15 2.86

20 3.44

25 3.87

30 4.15

35 4.37

40 4.51

45 4.58

50 4.02

55 4.64

江西理工大学上数值分析期末实验

hold on

plot(t,y2,'b')%绘制三次拟合函数图（其中y2是拟合之后的数据） xlabel('时间t（分钟）') %注释x轴 ylabel('含碳量/10^-4') %注释y轴 title('四次拟合图') %注释图名 grid%坐标系网格化

四、实验数据结果及分析

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

y1=0.00003436415436t3-0.00521556221556t2+0.26339852739853t+0.01783882783883

拟合函数与被拟合函数图之间的对比如下： （1） 红色星号为原始数据；

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

由此可见拟合函数与原函数离散数据点拟合成程度相当好，通过[p,s]=polyfit(t,y,n)对拟合误差进行分析，如图：

江西理工大学上数值分析期末实验

图2-2

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

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

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

y2=0.00000060256410t4-0.00003191789692 t3-0.00293227466977t2

+0.23806931494432t+0.06044871794872 拟合如图

2-3

图2-3

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

江西理工大学上数值分析期末实验

图2-4

由此可见，四次拟合误差0.49493<0.50824(三次拟合误差)，精度更高。

五、实验结论

在用高阶多项式对某一函数进行曲线拟合时, 并不是拟合出来的多项式与被拟合函数在整个区间上都能符合,polyfit()只能保证在输入数据x所能达到的区间上及其附近.求得的多项式可以最大限度在逼近原函数。利用最小二乘法对本问题进行的曲线拟合精度较高，而且，在一般情况下，拟合的多项式次数越多，精度越高。

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

一、问题提出

计算下列积分值：

（1

）I=

01

(I 1.5343916)

（2）I=

sinx

(f(0) 1,I 0.9460831) x0

1

ex

（3） I= 2

4 x0

江西理工大学上数值分析期末实验

（4） I=

ln(1 x)

2 1 x0

1

要求：

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

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

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

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

二、问题分析

由上可知这四个积分找不到用初等函数表示的原函数，直接计算起来很困难，因此我们考虑利用函数在若干点得函数值，近似地计算该函数在一个区间上得定积分。这里采用的方法有三种：复合梯形公式，复合Simpson公式，Romberg算法。

三、实验程序及注释

1、复合梯形公式MATLAB程序： function I=T_quad(x,y)

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

% x为向量，被积函数自变量的等距节点； % y为向量，被积函数节点处的函数值； n=length(x);m=length(y); if n~=m

error('the length of X and Y must be equal!'); return; end

h=(x(n)-x(1))/(n-1); a=[1 2*ones(1,n-2) 1]; I=h/2 * sum(a.*y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2、复合Simpson公式MATLAB程序： function I=S_quad(x,y)

% x为向量，被积函数自变量的等距节点； % y为向量，被积函数节点处的函数值； n=length(x);m=length(y); if n~=m

error('the length of X and Y must be equal!'); return; end

if rem(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);

江西理工大学上数值分析期末实验

for k=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

I=h/6*sum(a.*y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3、Romberg算法MATLAB程序： function I=R_quad_iter(fun,a,b,ep) % Romberg求积公式，其中， % fun为被积函数；

% a,b为积分区间端点，要求a<b; % ep精度系数，缺省值为1e-5. if nargin < 4 ep=1e-5; end

m=1;h=b-a;

I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I; while 1

N=2^(m-1);h=h/2;I=I/2; for i=1:N;

I=I+h*feval(fun,a+(2*i-1)*h); end

T(m+1,1)=I;M=2*N;k=1; while M>1;

T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k)/(4^k-1)); M=M/2;k=k+1; end

if abs(T(k,k)-T(k-1,k-1))<ep break; end

m=m+1; end

I=T(k,k);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4、自适应步长梯形公式：

function I=R_quad_iter(fun,a,b,ep) % 梯形递推求积公式，其中， % fun为被积函数；

% a,b为积分区间端点，要求a<b; % ep精度系数，缺省值为1e-5. if nargin < 4 ep=1e-5;end; N=1;h=b-a;

T=h/2*(feval(fun,a)+feval(fun,b)); while 1

江西理工大学上数值分析期末实验

h=h/2;I=T/2; for k=1:N;

I=I+h*feval(fun,a+(2*k-1)*h); end

if abs(I-T)<ep break; end

N=2*N;T=I; end

对四个积分求解：

%

对I=

(I 1.5343916)求解

x=0:1/40:1/4;

y=sqrt(4-(sin(x)).^2); format long

I=T_quad(x,y) %调用复化梯形公式求得积分值I =0.49870482652029 I=S_quad(x,y) %调用复化辛普森公式求得积分值I =0.49871111844568 fun=inline('sqrt(4-(sin(x)).^2)');

I=R_quad_iter(fun,0,1) %调用龙贝格算法可得积分值I=0.498711118445678 x=0:1/400:1/4;

y=sqrt(4-(sin(x)).^2);

I=T_quad(x,y) %调用复化梯形公式求得积分值I = 0.49871105466684 I=S_quad(x,y) %调用复化辛普森公式求得积分值I = 0.49871111757532

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对I=

sinx

(f(0) 1,I 0.9460831)求解： x0

1

x=0.1:0.1:1;

y=sin(x)./x %可以得到y值，由题可知f(0)=1。

y=[1 0.99833416646828 0.99334665397531 0.98506735553780 0.97354585577163 0.95885107720841 0.94107078899173 0.92031098176813 0.89669511362440 0.87036323291943 0.84147098480790]; x=0:0.1:1;

I=T_quad(x,y) %调用复化梯形公式可以求得I=0.94583207186691

I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.94608316883807 x=0.01:0.01:1;%x使用不同的步长 y=sin(x)./x

I=T_quad(x,y) %调用复化梯形公式可以求得I=0.94608266887360

I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.94608313833507

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

ex

%对I= 求解： 2

4 x0

1

江西理工大学上数值分析期末实验

x=0:0.1:1;

y=exp(x)./(4+x.^2);

I=T_quad(x,y) %调用复化梯形公式可求得积分值I =0.39087531274504 I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.39081195686445 x=0:0.01:1;%x使用不同的步长 y=exp(x)./(4+x.^2);

I=T_quad(x,y) %调用复化梯形公式可求得积分值I= 0.39081184557538 I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.39081184557538

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对I=

ln(1 x)

求积分： 2 1 x0

1

x=0:0.1:1;

y=log(1+x)./(1+x.^2);

I=T_quad(x,y) %调用复化梯形公式可求得积分值I =0.27128371750865 I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.27220124573417 x=0:0.01:1;%x使用不同的步长 y=log(1+x)./(1+x.^2);

I=T_quad(x,y) %调用复化梯形公式可求得积分值I= 0.27218912310178 I=S_quad(x,y) %调用复化辛普森公式求得积分值I =0.27219826157968

%给定精度要求ε，试用变步长算法，确定最佳步长： fun=inline('sqrt(4-(sin(x)).^2)'); I=R_quad_iter(fun,0,1/4)

二、应用题

1.文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（一天文单位为地球到太阳的平均距离：9300万里）

在五个不同的时间对小行星作了五次观察，测得轨道上五个点的坐标数据如下表所示：

由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：

现需要建立椭圆的方程以供研究。

江西理工大学上数值分析期末实验

（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写

出线性方程组AX = b。

（2）用MATLAB求低价方程组的指令A \ b求出待定系数

。

（3）卫星轨道是一个椭圆，其周长的计算公式是:

c

s 4a1 sin2 d

a

2

式中，a是椭圆的半长轴，

距离，

径。

是地球中心与轨道中心（椭圆中心）的

。其中h为近地点距离，H为远地点距离，R = 6371(km)为地球半

有一颗人造卫星近地点距离h = 439 (km)，远地点距离H = 2384(km)。试分别按下列方案计算卫星轨道的周长，误差限取为

。

三、要求

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

2、 对基本题，分别取不同步长h (b a)/n，试比较计算结果（如n = 10, 20

等）， 并比较其结果；

4、 对应用题，用给定精度ε，试用(1)用逐次分半梯形法。(2)用逐次分半辛普生法，并确定最佳步长。

逐次梯形法的主函数程序：

function [jifen,num]=zhuci_tixing(a,b,tol) if (nargin==3) eps=1.0e-3; end delt=1; n=1; h=b-a;

T=h*(zhuci_tixing_f(a)+zhuci_tixing_f(b))/2; while delt>3*tol h=h/2;

江西理工大学上数值分析期末实验

T0=T; s=0; for i=1:n

x=a+h*(2*i-1); s=s+zhuci_tixing_f(x); end

T=T0/2+h*s; n=2*n; delt=T-T0; end jifen=T;

num=n; end

执行程序后的结果：

江西理工大学上数值分析期末实验

实验分析：

逐次分半的积分算法具有很高的精度，对于复杂的工程实践问题具有很高的应用价值。 同时，在利用复合梯形公式、复合Simpson公式以及Romberg算法等计算数值积分时，精度已经很高，并且随着步长越小，积分精度越高。另外，当给定一个精度的时候，我们还可以利用自适应步长法，确定所需要的最佳步长和结果。

同样在对比中可见：simpson公式的收敛速度比梯形公式的收敛速度快。