实验报告—数值分析

《数值分析》实验报告姓名：学号：专业：指导教师：刘建生教授日期：2015年12月25日实验一Lagrange/newton插值一：对于给定的一元函数)(xfy的n+1个节点值(),0,1,,jjyfxjn。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：求五次L计算(0.596)f,(0.99)f的值（提示：结果为(0.596)0.625732f,(0.99)1.05423f）jx1234567jy0.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式6L()x，计算的(1.8)f,(6.15)f值。（提示：结果为(1.8)0.164762f,(6.15)0.001266f）二：实验程序及注释MATLAB程序：functionf=lagrange(x0,y0,x)n=length(x0);m=length(y0);formatlongs=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(x-x0(j))/(x0(k)-x0(j));endends=s+y0(k)*p;Endf=s;endjx0.40.550.650.800.951.05jy0.410750.578150.696750.901.001.25382结果运行：结果与提示值完全吻合，说明Lagrange插值多项式的精度是很高的；)45)(35)(25)(15)(05()4)(3)(2)(1)(0()50)(40)(30)(20)(10()5)(4)(3)(2)(1()(fxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx同时，若采用三点插值和两点插值的方法，用三点插值的精度更高。若同时采用两点插值，选取的节点距离x越近，精度越高。三：采用newton插值进行计算算法程序如下：formatlong;x0=[0.40.550.650.800.951.05];y0=[0.410750.578150.696750.901.001.25382];x=0.596;n=max(size(x0));y=y0(1);%disp(y);s=1;dx=y0;fori=1:n-1dx0=dx;forj=1:n-idx(j)=(dx0(j+1)-dx0(j))/(x0(i+j)-x0(j));enddf=dx(1);s=s*(x-x0(i));y=y+s*df;%计算%%disp(y);enddisp(y)运行结果：绘制出曲线图：与结果相吻合。所以newton法和Lagrange法的思想是一样的。Lagrange适合理论分析，但Lagrange法不如newton法灵活。Lagrange如果节点个数改变，算法需要重新编写，而Newton法克服这一缺点，所以应用更为灵活。实验二函数逼近与曲线拟合一、问题提出在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。要求：1、用最小二乘法进行曲线拟合；2、近似解析表达式为f(x)=a1t+a2t2+a3t3；3、计算出拟合函数f(x)，并列出出f(x)与y(x)的误差；4、另外选取一个近似表达式，尝试拟合效果的比较；5、绘制出曲线拟合图。二、问题分析三、从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。三、实验程序及注释三次拟合程序（最小二乘法）：t=[0510152025303540455055]%输入时间t的数据y=[01.272.162.863.443.874.154.374.514.584.024.64]%输入含碳量数据[p,s]=polyfit(t,y,3)%调用MATLAB最小二乘法的程序进行三次拟合并给出误差分析formatlong%14位精度小数plot(t,y,'*r')%绘制被拟合数据点的离散图holdonplot(t,y1,'b')%绘制三次拟合函数图（其中y1是拟合之后的数据）xlabel('时间t（分钟）')%注释x轴ylabel('含碳量/10^-4')%注释y轴title('三次拟合图')%注释图名grid%坐标系网格化四次拟合程序（最小二乘法）：[p,s]=polyfit(t,y,4)%调用MATLAB最小二乘法的程序进行四次拟合并给出误差分析formatlong%14位精度小数plot(t,y,'*r')%绘制被拟合数据点的离散图t(分钟)0510152025303540455055y（104）01.272.162.863.443.874.154.374.514.584.024.64holdonplot(t,y2,'b')%绘制三次拟合函数图（其中y2是拟合之后的数据）xlabel('时间t（分钟）')%注释x轴ylabel('含碳量/10^-4')%注释y轴title('四次拟合图')%注释图名grid%坐标系网格化四、实验数据结果及分析三次拟合可以得到其拟合多项式为：1y=0.00003436415436t3-0.00521556221556t2+0.26339852739853t+0.01783882783883拟合函数与被拟合函数图之间的对比如下：（1）红色星号为原始数据；（2）带圈的曲线为最小二乘后而成的结果曲线。由此可见拟合函数与原函数离散数据点拟合成程度相当好，通过[p,s]=polyfit(t,y,n)对拟合误差进行分析，如图：图2-2由此可知，三次拟合精度较好。为了提高结果的可信度，我们另外选取一个近似表达式，尝试拟合效果的比较。于是，进行四次拟合：其中，拟合得到的多项式为：2y=0.00000060256410t4-0.00003191789692t3-0.00293227466977t2+0.23806931494432t+0.06044871794872拟合如图2-3图2-3同样对四次拟合进行误差分析可得：图2-4由此可见，四次拟合误差0.494930.50824(三次拟合误差)，精度更高。五、实验结论在用高阶多项式对某一函数进行曲线拟合时,并不是拟合出来的多项式与被拟合函数在整个区间上都能符合,polyfit()只能保证在输入数据x所能达到的区间上及其附近.求得的多项式可以最大限度在逼近原函数。利用最小二乘法对本问题进行的曲线拟合精度较高，而且，在一般情况下，拟合的多项式次数越多，精度越高。实验三数值积分与数值微分一、问题提出计算下列积分值：（1）1/420I=4-sin(1.5343916)xdxI（2）10sinI=((0)1,0.9460831)xdxfIx（3）120I=4xedxx（4）120ln(1)I=1xdxx要求：1、编制数值积分算法的程序；2、分别用两种算法计算同一个积分，并比较其结果；3、分别取不同步长()/hban，试比较计算结果（如n=10,20等）；4、给定精度要求ε，试用变步长算法，确定最佳步长。二、问题分析由上可知这四个积分找不到用初等函数表示的原函数，直接计算起来很困难，因此我们考虑利用函数在若干点得函数值，近似地计算该函数在一个区间上得定积分。这里采用的方法有三种：复合梯形公式，复合Simpson公式，Romberg算法。三、实验程序及注释1、复合梯形公式MATLAB程序：functionI=T_quad(x,y)%复化梯形求积公式，其中，%x为向量，被积函数自变量的等距节点；%y为向量，被积函数节点处的函数值；n=length(x);m=length(y);ifn~=merror('thelengthofXandYmustbeequal!');return;endh=(x(n)-x(1))/(n-1);a=[12*ones(1,n-2)1];I=h/2*sum(a.*y);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2、复合Simpson公式MATLAB程序：functionI=S_quad(x,y)%x为向量，被积函数自变量的等距节点；%y为向量，被积函数节点处的函数值；n=length(x);m=length(y);ifn~=merror('thelengthofXandYmustbeequal!');return;endifrem(n-1,2)~=0%如果n-1不能被2整除，则调用复化梯形公式I=T_quad(x,y);return;endN=(n-1)/2;h=(x(n)-x(1))/N;a=zeros(1,n);fork=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4;a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3、Romberg算法MATLAB程序：functionI=R_quad_iter(fun,a,b,ep)%Romberg求积公式，其中，%fun为被积函数；%a,b为积分区间端点，要求ab;%ep精度系数，缺省值为1e-5.ifnargin4ep=1e-5;endm=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while1N=2^(m-1);h=h/2;I=I/2;fori=1:N;I=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;whileM1;T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k)/(4^k-1));M=M/2;k=k+1;endifabs(T(k,k)-T(k-1,k-1))epbreak;endm=m+1;endI=T(k,k);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%4、自适应步长梯形公式：functionI=R_quad_iter(fun,a,b,ep)%梯形递推求积公式，其中，%fun为被积函数；%a,b为积分区间端点，要求ab;%ep精度系数，缺省值为1e-5.ifnargin4ep=1e-5;end;N=1;h=b-a;T=h/2*(feval(fun,a)+feval(fun,b));while1h=h/2;I=T/2;fork=1:N;I=I+h*feval(fun,a+(2*k-1)*h);endifabs(I-T)epbreak;endN=2*N;T=I;end对四个积分求解：%对1/420I=4-sin(1.5343916)xdxI求解x=0:1/40:1/4;y=sqrt(4-(sin(x)).^2);formatlongI=T_quad(x,y)%调用复化梯形公式求得积分值I=0.49870482652029I=S_quad(x,y)%调用复化辛普森公式求得积分值I=0.49871111844568fun=inline('sqrt(4-(sin(x)).^2)');I=R_quad_iter(fun,0,1)%调用龙贝格算法可得积分值I=0.498711118445678x=0:1/400:1/4;y=sqrt(4-(sin(x)).^2);I=T_quad(x,y)%调用复化梯形公式求得积分值I=0.49871105466684I=S_quad(x,y)%调用复化辛普森公式求得积分值I=0.49871111757532%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对10sinI=((0)1,0.9460831)xdxfIx求解：x=0.1:0.1:1;y=sin(x)./x%可以得到y值，由题可知f(0)=1。y=[10.998334166468280.993346653975310.985067355537800.973545855771630.958851077208410.941070788991730