实验报告数值分析.doc

数值分析实验报告姓 名： 学 号： 专 业： 指导教师： 实验一 Lagrange/newton插值一：对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下： 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.578150.696750.90 1.00 1.25382 求五次L计算,的值（提示：结果为, ） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，计算的,值。（提示：结果为, ） 二：实验程序及注释MATLAB程序：function f=lagrange(x0,y0,x )n=length(x0);m=length(y0); format longs=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插值多项式的精度是很高的；同时，若采用三点插值和两点插值的方法，用三点插值的精度更高。若同时采用两点插值，选取的节点距离x越近，精度越高。三：采用newton插值进行计算算法程序如下：format long; x0=0.4 0.55 0.65 0.80 0.95 1.05 ;y0=0.41075 0.578150.696750.90 1.00 1.25382 ;x=0.596; 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);enddisp(y)运行结果：绘制出曲线图：与结果相吻合。所以newton法和Lagrange法的思想是一样的。Lagrange适合理论分析，但Lagrange法不如newton法灵活。Lagrange如果节点个数改变，算法需要重新编写，而Newton法克服这一缺点，所以应用更为灵活。实验二 函数逼近与曲线拟合 一、问题提出 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。 t(分钟)0510152025303540455055y（10）01.272.162.863.443.874.154.374.514.584.024.64要求： 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为f(x)=at+ at+ at；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 onplot(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')%绘制被拟合数据点的离散图hold onplot(t,y2,'b')%绘制三次拟合函数图（其中y2是拟合之后的数据）xlabel('时间t（分钟）') %注释x轴ylabel('含碳量/10-4') %注释y轴title('四次拟合图') %注释图名grid%坐标系网格化四、实验数据结果及分析三次拟合可以得到其拟合多项式为：=0.00003436415436t-0.00521556221556t+0.26339852739853t+0.01783882783883 拟合函数与被拟合函数图之间的对比如下：（1） 红色星号为原始数据；（2） 带圈的曲线为最小二乘后而成的结果曲线。由此可见拟合函数与原函数离散数据点拟合成程度相当好，通过p,s=polyfit(t,y,n)对拟合误差进行分析，如图：图2-2由此可知，三次拟合精度较好。为了提高结果的可信度，我们另外选取一个近似表达式，尝试拟合效果的比较。于是，进行四次拟合：其中，拟合得到的多项式为：=0.00000060256410t-0.00003191789692 t-0.00293227466977t+0.23806931494432t+0.06044871794872拟合如图2-3图2-3同样对四次拟合进行误差分析可得：图2-4由此可见，四次拟合误差0.49493<0.50824(三次拟合误差)，精度更高。五、实验结论在用高阶多项式对某一函数进行曲线拟合时,并不是拟合出来的多项式与被拟合函数在整个区间上都能符合,polyfit()只能保证在输入数据x所能达到的区间上及其附近.求得的多项式可以最大限度在逼近原函数。利用最小二乘法对本问题进行的曲线拟合精度较高，而且，在一般情况下，拟合的多项式次数越多，精度越高。实验三 数值积分与数值微分一、问题提出 计算下列积分值：（1） （2） （3） （4） 要求：1、 编制数值积分算法的程序； 2、 分别用两种算法计算同一个积分，并比较其结果； 3、 分别取不同步长，试比较计算结果（如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;endh=(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;endif rem(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);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;endI=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;endm=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)=(4k*T(m+1,k)-T(m,k)/(4k-1); M=M/2;k=k+1; end if abs(T(k,k)-T(k-1,k-1)<ep break; end m=m+1;endI=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对四个积分求解：%对求解x=0:1/40:1/4;y=sqrt(4-(sin(x).2);format longI=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 %对求解： 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.94583207186691I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.94608316883807x=0.01:0.01:1;%x使用不同的步长y=sin(x)./xI=T_quad(x,y) %调用复化梯形公式可以求得I=0.94608266887360I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.94608313833507%对求解：x=0:0.1:1;y=exp(x)./(4+x.2);I=T_quad(x,y) %调用复化梯形公式可求得积分值I =0.39087531274504I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.39081195686445x=0:0.01:1;%x使用不同的步长y=exp(x)./(4+x.2);I=T_quad(x,y) %调用复化梯形公式可求得积分值I= 0.39081184557538I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.39081184557538%对求积分：x=0:0.1:1;y=log(1+x)./(1+x.2);I=T_quad(x,y) %调用复化梯形公式可求得积分值I =0.27128371750865I=S_quad(x,y) %调用复化辛普森公式求得积分值I= 0.27220124573417x=0:0.01:1;%x使用不同的步长y=log(1+x)./(1+x.2);I=T_quad(x,y) %调用复化梯形公式可求得积分值I= 0.27218912310178I=S_quad(x,y) %调用复化辛普森公式求得积分值I =0.27219826157968%给定精度要求，试用变步长算法，确定最佳步长：fun=inline('sqrt(4-(sin(x).2)');I=R_quad_iter(fun,0,1/4)二、应用题1.文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（一天文单位为地球到太阳的平均距离：9300万里）在五个不同的时间对小行星作了五次观察，测得轨道上五个点的坐标数据如下表所示：   P1P2P3P4P5x坐标5.7646.2866.7597.1687.408y坐标0.6481.2021.8232.5267.408由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：现需要建立椭圆的方程以供研究。 （1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组AX = b。（2）用MATLAB求低价方程组的指令A b求出待定系数 。（3）卫星轨道是一个椭圆，其周长的计算公式是:式中，a是椭圆的半长轴， 是地球中心与轨道中心（椭圆中心）的距离， 。其中h为近地点距离，H为远地点距离，R = 6371(km)为地球半径。   有一颗人造卫星近地点距离h = 439 (km)，远地点距离H = 2384(km)。试分别按下列方案计算卫星轨道的周长，误差限取为 。三、要求 1、 编制数值积分算法的程序； 2、 对基本题，分别取不同步长，试比较计算结果（如n = 10, 20等）， 并比较其结果； 4、 对应用题，用给定精度，试用(1)用逐次分半梯形法。(2)用逐次分半辛普生法，并确定最佳步长。 逐次梯形法的主函数程序：function jifen,num=zhuci_tixing(a,b,tol)if (nargin=3) eps=1.0e-3;enddelt=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;endjifen=T;num=n;end执行程序后的结果：实验分析： 逐次分半的积分算法具有很高的精度，对于复杂的工程实践问题具有很高的应用价值。同时，在利用复合梯形公式、复合Simpson公式以及Romberg算法等计算数值积分时，精度已经很高，并且随着步长越小，积分精度越高。另外，当给定一个精度的时候，我们还可以利用自适应步长法，确定所需要的最佳步长和结果。 同样在对比中可见：simpson公式的收敛速度比梯形公式的收敛速度快。实验四 线性方程组的直接解法一、问题提出 给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、 设线性方程组 参考解：2、 设对称正定阵系数阵线方程组 参考解：3，三对角形线性方程组 参考解：要求：1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）； 2、 应用结构程序设计编出通用程序； 3、 比较计算结果，分析数值解误差的原因； 4、 尽可能利用相应模块输出系数矩阵的三角分解式。二、问题分析直接法就是经过有限步算术运算，可求得线性方程组精确解得方法（若计算过程中没有舍入误差）。但是实际运算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解。我们经常采用的一些方法，如Gauss顺序消去法与Gauss列主元消去法、平方根法与改进平方根法、追赶法等等。三、 实验程序及注释 列主元Gauss消去法：function x,det,index=Gauss(A,b)% 求线性方程的列主元Gauss消去法% A 为方程组的系数矩阵% b 为方程组的右端项% x 为方程组的解% det 为系数矩阵A的行列式的值% index 为指标变量，index=0 表示计算失败，index=1 表示计算成功% det 为系数矩阵 A 的行列式的值n,m=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出错误的信息if n=m error('A的行与列必须相等！'); return;end% 当方程组与右端项的维数不匹配时，停止计算，输出错误信息if m=nb error('the columns of A must be equal the length of b'); return;end% 开始计算，先赋初值index=1;det=1;x=zeros(n,1);for k=1:n-1 % 选主元 a_max=0; for i=k:n if abs(A(i,k)>a_max a_max=abs(A(i,k);r=i; end end if a_max<1e-10 index=0;return; end % 交换两行 if r>k; for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end % 消元过程 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end det=det*A(k,k);enddet=det*A(n,n);% 回代过程if abs(A(n,n)<1e-10 index=0;return;endfor k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k);end改进平方根法：function L,D,index=LDL_Decom(A)% 求正定对称矩阵A的改进平方根分解，也称LDLT分解% A为要分解的矩阵% L为下三角阵% D为对角阵% index 为指标变量，index=0 表示计算失败，index=1 表示计算成功n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);T=zeros(n);index=1;for k=1:n d(k)=A(k,k); for j=1:k-1 d(k)=d(k)-L(k,j)*T(k,j); end if abs(d(k)<1e-10 index=0;return; end for i=k+1:n T(i,k)=A(i,k); for j=1:k-1 T(i,k)=T(i,k)-T(i,j)*L(k,j); end L(i,k)=T(i,k)/d(k); endend D=diag(d);追赶法：function x=followup(A,b)% 采用追赶法求解% 系数矩阵：A% 常数向量：b% 线性方程组的解：xn=rank(A);for (i=1:n) if (A(i,i)=0) disp('error:对角元素为0'); return; endend; d=ones(n,1); a=ones(n-1,1); c=ones(n-1); for(i=1:n-1) a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i); end d(n,1)=A(n,n); for(i=2:n) d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(i-1,1); end x(n,1)=b(n,1)/d(n,1); for(i=(n-1):-1:1) x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1); end第一问：A=4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1b=5 12 3 2 3 46 13 38 19 -21'x,det,index=Gauss(A,b)L,U,P=lu(A)第二问：A=4 2 -4 0 2 4 0 0;2 2 -1 -2 1 3 2 0;-4 -1 14 1 -8 -3 5 6;0 -2 1 6 -1 -4 -3 3;2 1 -8 -1 22 4 -10 -3;4 3 -3 -4 4 11 1 -4;0 2 5 -3 -10 1 14 2;0 0 6 3 -3 -4 2 19b=0 -6 20 23 9 -22 -15 45'L,D,index=LDL_Decom(A)R=L'*bm=D*L'mrL,U,P=lu(A)第三问：A=4 -1 0 0 0 0 0 0 0 0;-1 4 -1 0 0 0 0 0 0 0;0 -1 4 -1 0 0 0 0 0 0;0 0 -1 4 -1 0 0 0 0 0 ;0 0 0 -1 4 -1 0 0 0 0;0 0 0 0 -1 4 -1 0 0 0;0 0 0 0 0 -1 4 -1 0 0;0 0 0 0 0 0 -1 4 -1 0;0 0 0 0 0 0 0 -1 4 -1;0 0 0 0 0 0 0 0 -1 4b=7 5 -13 2 6 -12 14 -4 5 -5'x=followup(A,b)L,U,P=lu(A)四、实验数据结果及分析1、第一问中的方程组属于一般的线性方程组，这里采用列主元高斯消去法，求得对比精确值，我们发现差距很大。估计是所给的参考值错误。同样利用lu()函数可以系数矩阵A的三角分解式：2、 该方程组为三对角形线性方程组，可以采用追赶法解决类似的问题。得到的结果对比题中的参考值，精度很高。三角分解为：五、实验结论本实验运用最基本的算法（高斯消去法以及某些变形）。这类方法是解低阶稠密矩阵方程组以及某些大型稀疏矩阵方程组的有效方法。列主元消去法比Gauss顺序消去法更加能保证结果精度，因为列主元消去法始终保持除数最大，避免了因为计算机内存结构而引起的“大数”吃“小数”，故在工程中比较常用。实验五 解线性方程组的迭代法一、问题提出 对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。实验程序代码：一，%jacobi迭代法计算线性方程组function x,k=jacobi_f(A,b,x0,tol)max1=300; D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); B=D(L+U); f=Db; x=B*x0+f; k=1; while norm(x-x0)>=tol x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代次数超过300次，方程组可能无解'); return; end二，%高斯赛德尔迭代法计算线性方程组function x,k=guass_seid_f(A,b,x0,tol) max1=300; D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); G=(D-L)U; f=(D-L)b; x=G*x0+f; k=1; while norm(x-x0)>=tol x0=x; x=G*x0+f; k=k+1; if(k>=max1) disp('迭代次数太多，可能不收敛') end %k,x' end三，%超松弛迭代法（sor)方法函数%超松弛迭代法的统一方程 function x,k=sor_f(A,b,x0,w,tol)max=300;if (w<=0|w>=2) error; return;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-w*L)*b;x=B*x0+f;k=1;while norm(x-x0)>=tol x0=x; x=B*x0+f; k=k+1; if (k>max) disp('迭代次数过多，sor方法可能不收敛'); return; end %k,x' %输出每一步迭代的结果，可不要End实验结果：对比题中的参考值，精度很高。实验六 非线性方程求根一、问题提出 设方程有三个实根现采用下面六种不同计算格式，求 f(x)=0的根或1、 2、 3、 4、 5、 6、 实验程序代码： 本程序采用牛顿法求解非线性方程组： 也可自行编写简单迭代算法进行求解，此处没有给出。程序代码如下：%newton迭代法求解非线性方程function gen,time =newton(f,x0,tol)if (nargin=2) tol=1.0e-5;end df=diff(f); x1=x0; time=0; wucha=0.1; while(wucha>tol) time=time+1; fx=subs(f,x1); df=subs(df,x1); gen=x1-fx/df; wucha=abs(gen-x1); x1=gen; endend分别带入各个迭代格式得：由收敛准则: 知道，当且仅当迭代方程满足以上的收敛条件时，非线性方程组才有解。实验七 矩阵特征值问题计算一、问题提出 利用冪法或反冪法，求方阵的按模最大或按模最小特征值及其对应的特征向量。 设矩阵A的特征分布为： 且 试求下列矩阵之一 （1） 求，及 取 结果（2） 求及 取 结果：（3） 求及取结果 (4) 取 这是一个收敛很慢的例子,迭代次才达到结果 （5） 有一个近似特征值,试用幂法求对应的特征向量，并改进特征值（原点平移法）。取 结果 运用MATLAB指令：V,D=eig(A)； 也可用幂法求解，程序在下面给出，可自行调试进行求解：（1）（2）（3）下面给出幂法的算法程序：已调试：function namda,q,time,date_na,date_q=mifa(A,tol,x0) if nargin =1 tol=1e-7; temp=length(A); x0=ones(temp,1); end if nargin=2 temp=length(A); x0=ones(temp,1); end m1=0; u=x0; time=0;while time<500 v=A*u; vmax,i=max(abs(v); m=v(i); u=v/m; if abs(m-m1)<tol break; end m1=m; time=time+1; date_na(time,1)=m; date_q(:,time)=u;endnamda=m;q=u;end 实验八 常微分方程初值问题数值解法1、问题的提出 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法(Euler 预测校正法)，四阶Rung-Kutta方法求其数值解，诸如以下问题： 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。 2、 问题的分析（1）、利用Euler公式：当取定初值y(0)=3时，即可由上式递推计算出y1,y2,.,y(n).（2） 、利用改进的Euler公式：（3） 利用四阶Rung-Kutta公式：3、实验主程序及部分注释 3.1 Euler法编程主程序:%Euler法 function x,y=Euler(fun,span,y0,h)x=span(1):h:span(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(fun,x(n),y(n);endR=x' y' 3.2 改进Euler法编程主程序:%改进的Euler法(Euler 预测校正法)function x,y = Euler1(fun,span,y0,h)x=span(1):h:span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); k2