《数值分析》课程实验报告.docx
- 文档编号:26968052
- 上传时间:2023-06-24
- 格式:DOCX
- 页数:12
- 大小:23.63KB
《数值分析》课程实验报告.docx
《《数值分析》课程实验报告.docx》由会员分享,可在线阅读,更多相关《《数值分析》课程实验报告.docx(12页珍藏版)》请在冰豆网上搜索。
《数值分析》课程实验报告
《数值分析》课程实验报告
《数值分析》课程实验报告姓名:
学号:
学院:
机电学院日期:
2015年X月X日目录实验一函数插值方法1实验二函数逼近与曲线拟合5实验三数值积分与数值微分7实验四线方程组的直接解法9实验五解线性方程组的迭代法15实验六非线性方程求根19实验七矩阵特征值问题计算21实验八常微分方程初值问题数值解法24实验一函数插值方法一、问题提出对于给定的一元函数的n+1个节点值。
试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。
数据如下:
(1)
0.40.550.650.800.951.050.410750.578150.696750.901.001.25382求五次Lagrange多项式,和分段三次插值多项式,计算,的值。
(提示:
结果为,)
(2)
12345670.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式,计算的,值。
(提示:
结果为,)
二、要求1、利用Lagrange插值公式编写出插值多项式程序;
2、给出插值多项式或分段三次插值多项式的表达式;
3、根据节点选取原则,对问题
(2)用三点插值或二点插值,其结果如何;
4、对此插值问题用Newton插值多项式其结果如何。
Newton插值多项式如下:
其中:
三、目的和意义1、学会常用的插值方法,求函数的近似表达式,以解决其它实际问题;
2、明确插值多项式和分段插值多项式各自的优缺点;
3、熟悉插值方法的程序编制;
4、如果绘出插值函数的曲线,观察其光滑性。
四、实验步骤
(1)
0.40.550.650.800.951.050.410750.578150.696750.901.001.25382求五次Lagrange多项式,和分段三次插值多项式,计算,的值。
(提示:
结果为,)
第一步:
先在matlab中定义lagran的M文件为拉格朗日函数代码为:
function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:
n+1v=1;forj=1:
n+1if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:
)=v;endc=y*l;end第二步:
然后在matlab命令窗口输入:
x=[0.40.550.650.80,0.951.05];y=[0.410750.578150.696750.901.001.25382];lagran(x,y)回车得到:
ans=121.6264-422.7503572.5667-377.2549121.9718-15.0845由此得出所求拉格朗日多项式为p(x)=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845第三步:
在编辑窗口输入如下命令:
x=[0.40.550.650.80,0.951.05];y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845;plot(x,y)命令执行后得到如下图所示图形,然后x=0.596;y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084y=0.6262得到f(0.596)=0.6262同理得到f(0.99)=1.0547
(2)
12345670.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式,和分段三次插值多项式,计算的,值。
(提示:
结果为,)
实验步骤:
第一步定义function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:
n+1v=1;forj=1:
n+1if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:
)=v;endc=y*l;end定义完拉格朗日M文件第二步:
然后在matlab命令窗口输入:
x=[1234567];y=[0.3680.1350.0500.0180.0070.0020.001];lagran(x,y)回车得到:
ans=0.0001-0.00160.0186-0.11750.4419-0.96830.9950由此得出所求拉格朗日多项式为p(x)=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950第三步:
在编辑窗口输入如下命令:
x=[1234567];y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950;plot(x,y)命令执行后得到如下图所示图形,然后x=1.8;y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084y=0.1650得到f(0.596)=0.6262同理得到f(6.15)=2.3644五、实验结论插值是在离散数据的基础上补插连续函数,使得这条连续曲线通过全部给定的离散数据点,它是离散函数逼近的重要方法,利用它可通过函数在有限个点处的取值状况,估算出函数在其他点处的近似值。
实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量与时间t的拟合曲线。
t(分)051015202530354045505501.272.162.863.443.874.154.374.514.584.024.64二、要求1、用最小二乘法进行曲线拟合;
2、近似解析表达式为;
3、打印出拟合函数,并打印出与的误差,;
4、另外选取一个近似表达式,尝试拟合效果的比较;
5、*绘制出曲线拟合图。
三、目的和意义1、掌握曲线拟合的最小二乘法;
2、最小二乘法亦可用于解超定线代数方程组;
3、探索拟合函数的选择与拟合精度间的关系四、实验步骤:
第一步先写出线性最小二乘法的M文件functionc=lspoly(x,y,m)n=length(x);b=zeros(1:
m+1);f=zeros(n,m+1);fork=1:
m+1f(:
k)=x.^(k-1);enda=f'*f;b=f'*y';c=a\b;c=flipud(c);第二步在命令窗口输入:
lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)回车得到:
ans=-0.00240.20370.2305即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305在编辑窗口输入如下命令:
x=[0,5,10,15,20,25,30,35,40,45,50,55];y=-0.0024*x.^2+0.2037*x+0.2305;plot(x,y)命令执行得到如下图五、实验结论分析复杂实验数据时,常采用分段曲线拟合方法。
利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。
分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。
实验三数值积分与数值微分一、问题提出选用复合梯形公式,复合Simpson公式,Romberg算法,计算
(1)
(2)
(3)
(4)
二、要求1、编制数值积分算法的程序;
2、分别用两种算法计算同一个积分,并比较其结果;
3、分别取不同步长,试比较计算结果(如n=10,20等);
4、给定精度要求ε,试用变步长算法,确定最佳步长。
三、目的和意义1、深刻认识数值积分法的意义;
2、明确数值积分精度与步长的关系;
3、根据定积分的计算方法,可以考虑二重积分的计算问题。
四、实验步骤第一步:
编写各种积分的程序复合梯形程序如下:
functionI=TX(x,y)n=length(x);m=length(y);ifn~=merror('ThelengthsofXandYmustbeequal');return;endh=(x(n)-x
(1))/(n-1);a=[12*ones(1,n-2)1];I=h/2*sum(a.*y);复合Simpson程序如下:
functions=simpr1(f,a,b,n)h=(b-a)/(2*n);s1=0;s2=0;fork=1:
10x=a+h*(2*k-1);s1=s1+feval(f,x);endfork=1:
(10-1)x=a+h*2*k;s2=s2+feval(f,x);ends=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;endRomberg程序如下:
functionI=Romber_yang(fun,a,b,ep)ifnargin4ep=1e-5;end;m=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:
NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;whileMT(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);第二步:
用复合梯形公式,复合Simpson公式,Romberg公式对各个积分进行计算1、对于积分Ι=0144-sin2xdx,梯形积分T=0.49871011,辛普森积分S=0.49871111,Romberg积分R=0.49871111。
2、对于积分Ι=01sinXXdx,f(0)=1,梯形积分T=0.94607307,辛普森积分S=0.94607308,Romberg积分R=0.94607307。
3、对于积分Ι=01eX4+X2dx,梯形积分T=0.39081248,辛普森积分S=0.39081185,Romberg积分R=0.39081885。
4、对于积分Ι=01ln1+X1+X2dx,梯形积分T=0.27218912,辛普森积分S=0.27219844,Romberg积分R=0.27219827。
五、实验结论,通过本实验学会复合梯形公式,复合Simpson公式,Romberg公式的编程与应用,掌握MATLAB提供的计算积分的各种函数的使用方法。
实验四线方程组的直接解法一、问题提出给出下列几个不同类型的线性方程组,请用适当算法计算其解。
1、设线性方程组2、设对称正定阵系数阵线方程组3、三对角形线性方程组二、要求1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法;
平方根法与改进平方根法;
追赶法求解(选择其一);
2、应用结构程序设计编出通用程序;
3、比较计算结果,分析数值解误差的原因;
4、尽可能利用相应模块输出系数矩阵的三角分解式。
三、目的和意义1、通过该课题的实验,体会模块化结构程序设计方法的优点;
2、运用所学的计算方法,解决各类线性方程组的直接算法;
3、提高分析和解决问题的能力,做到学以致用;
4、通过三对角形线性方程组的解法,体会稀疏线性方程组解法的特点。
四、实验步骤:
列主元高斯消去法的matlab的M文件程序function[x,det,index]=Gauss(A,b)%求线形方程组的列主元Gauss消去法,其中,%A为方程组的系数矩阵;
%b为方程组的右端项;
%x为方程组的解;
%det为系数矩阵A的行列式的值;
%index为指标变量,index=0表示计算失败,index=1表示计算成功。
[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时,停止计算,并输出出错信息。
ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!
');return;end%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息ifm~=nberror('ThecolumnsofAmustbeequalthelengthofb!
');return;end%开始计算,先赋初值index=1;det=1;x=zeros(n,1);fork=1:
n-1%选主元a_max=0;fori=k:
nifabs(A(i,k))a_maxa_max=abs(A(i,k));r=i;endendifa_max1e-10index=0;return;end%交换两行ifrkforj=k:
nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;end%消元过程fori=k+1:
nm=A(i,k)/A(k,k);forj=k+1:
nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);%回代过程ifabs(A(n,n))1e-10index=0;return;endfork=n:
-1:
1forj=k+1:
nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end然后在命令窗口输入A=[42-3-1210000;86-5-3650100;42-2-132-1031;0-215-13-1194;-426-167-3323;86-8571726-35;02-13-425301;1610-11-917342-122;462-713920124;00-18-3-24-863-1];b=[51232346133819-21];gauss(A,b)ans=1.0000-1.00000.00001.00002.00000.00003.00001.0000-1.00002.0000高斯-约当消去法maltab的M文件程序function[x,flag]=Gau_Jor(A,b)%求线形方程组的列主元Gauss-约当法消去法,其中,%A为方程组的系数矩阵;
%b为方程组的右端项;
%x为方程组的解;
[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时,停止计算,并输出出错信息。
ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!
');return;end%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息ifm~=nberror('ThecolumnsofAmustbeequalthelengthofb!
');return;end%开始计算,先赋初值flag='ok';x=zeros(n,1);fork=1:
n%选主元max1=0;fori=k:
nifabs(A(i,k))max1max1=abs(A(i,k));r=i;endendifmax11e-10falg='failure';return;end%交换两行ifrkforj=k:
nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元过程b(k)=b(k)/A(k,k);forj=k+1:
nA(k,j)=A(k,j)/A(k,k);endfori=1:
nifi~=kforj=k+1:
nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend%输出xfori=1:
nx(i)=b(i);end然后保存后在命令窗口输入:
A=[42-402400;22-1-21320;-4-1141-8-356;0-216-1-4-33;21-8-1224-10-3;43-3-44111-4;025-3-101142;0063-3-4219];b=[0-620239-22-1545];Gau_Jor(A,b)ans=121.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.0185五、实验结论用LU法,调用matlab中的函数lu中,L往往不是一个下三角,但可以直接计算不用它的结果来计算,不用进行行变换。
如果进行行变b也要变,这样会很麻烦。
实验五解线性方程组的迭代法一、问题提出对实验四所列目的和意义的线性方程组,试分别选用Jacobi迭代法,Gauss-Seidel迭代法和SOR方法计算其解。
二、要求1、体会迭代法求解线性方程组,并能与消去法做以比较;
2、分别对不同精度要求,如由迭代次数体会该迭代法的收敛快慢;
3、对方程组2,3使用SOR方法时,选取松弛因子ω=0.8,0.9,1,1.1,1.2等,试看对算法收敛性的影响,并能找出你所选用的松弛因子的最佳者;
4、给出各种算法的设计程序和计算结果。
三、目的和意义1、通过上机计算体会迭代法求解线性方程组的特点,并能和消去法比较;
2、运用所学的迭代法算法,解决各类线性方程组,编出算法程序;
3、体会上机计算时,终止步骤或k(给予的迭代次数),对迭代法敛散性的意义;
4、体会初始解,松弛因子的选取,对计算结果的影响。
四、实验步骤第一步编写实验所需的Jacobi迭代法,Gauss-Seidel迭代法,SOR迭代法的程序。
Jacobi迭代法:
function[x,k,index]=J(A,b,ep,itmax)ifnargin4itmax=100;endifnargin3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1fori=1:
ny(i)=b(i);forj=1:
nifj~=iy(i)=y(i)-A(i,j)*x(j);endendifabs(A(i,i))1e-10|k==itmaxindex=0;return;endy(i)=y(i)/A(i,i);endifnorm(y-x,inf)epbreak;endx=y;k=k+1;endGauss-Seidel迭代法:
function[x,k,index]=G(A,b,ep,itmax)ifnargin4itmax=100;endifnargin3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1y=x;fori=1:
nz=b(i);forj=1:
nifj~=iz=z-A(i,j)*x(j);endendifabs(A(i,i))1e-10|k==itmaxindex=0;return;endz=z/A(i,i);x(i)=z;endifnorm(y-x,inf)epbreak;endk=k+1;endSOR迭代法:
function[x,k,index]=SOR(A,b,ep,w,itmax)ifnargin5itmax=100;endifnargin4w=1;endifnargin3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1y=x;fori=1:
nz=b(i);forj=1:
nifj~=iz=z-A(i,j)*x(j);endendifabs(A(i,i))1e-10|k==itmaxindex=0;return;endz=z/A(i,i);x(i)=(1-w)*x(i)+w*z;endifnorm(y-x,inf)epbreak;endk=k+1;end第二步:
把实验给出的方程代入程序计算其结果。
1、设线性方程组2、设对称正定阵系数阵线方程组3、三对角形线性方程组五、实验结论迭代法是解线性方程组的一个重要的实用方法,特别适用于求解在实际中大量出现的,系数矩阵为稀疏阵的大型线性方程组。
通过此次实验学会了Jacobi迭代法,Gauss-Seidel迭代法,SOR迭代法的程序编写,并掌握了它们各自的优缺点及其适用条件。
实验六非线性方程求根一、问题提出设方程有三个实根现采用下面六种不同计算格式,求f(x)=0的根或1、2、3、4、5、6、二、要求1、编制一个程序进行运算,最后打印出每种迭代格式的敛散情况;
2、用事后误差估计来控制迭代次数,并且打印出迭代的次数;
3、初始值的选取对迭代收敛有何影响;
4、分析迭代收敛和发散的原因。
三、目的和意义1、通过实验进一步了解方程求根的算法;
2、认识选择计算格式的重要性;
3、掌握迭代算法和精度控制;
4、明确迭代收敛性与初值选取的关系。
四、实验步骤第一步:
编写实验所需的程序。
function[x_star,index,it]=DD(fun,x,ep,itmax)ifnargin4itmax=100;endifnargin3ep=1e-5;endindex=0;k=1;whilekitmaxx1=x;x=feval(fun,x);ifabs(x-x1)epindex=1;break;endk=k+1;endx_star=x;it=k;第二步:
分别用上述六种形式的表达式计算方程的根,结果如下。
1、,x1=0,x2=0。
2、,x1=无穷大,x2=-0.3473。
3、,x1=1.8794,x2=1.8794。
4、,x1=-0.3473,x2=-0.3473.。
5、,x1=1.8794,x2=1.8794。
6、,x1=1.8794,x2=-0.3473。
五、实验结论对于非线性方程,求它的解析解有时候是很困难的,但采用数值方法可以很容易地求它的近似解。
此次实验就是采用迭代法求非线性方程的根。
对于一个非线性方程,选用不同的迭代形式,因为其收敛程度不一样,造成其效率与精确度有很大的差别。
实验七矩阵特征值问题计算一、问题提出利用冪法或反冪法,求方阵的按模最大或按模最小特征值及其对应的特征向量。
设矩阵A的特征分布为:
且试求下列矩阵之一
(1)
求,及取结果
(2)
求及取结果:
(3)
求及取结果(4)取这是一个收敛很慢的例子,迭代次才达到结果(5)
有一个近似特征值,试用幂法求对应的特征向量,并改进特征值(原点平移法)。
取结果二、要求1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计;
2、会用原点平移法改进算法,加速收敛;
对矩阵B=A-PI取不同的P值,试求其效果;
3、试取不同的初始向量,观察对结果的影响;
4、对矩阵特征值的其它分布,如且如何计算。
三、目的和意义1、求矩阵的部分特征值问题具有重要实际意义,如求矩阵谱半径,稳定性问题往往归于求矩阵按模最小特征值;
2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧;
3、问题中的题(5),反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征向量。
四、实验步骤第一步:
写出实验所需的幂法求最大特征值及反幂法求
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析 数值 分析 课程 实验 报告