版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
数值分析课程实验指导书实验一函数插值方法一、问题提出对于给定的一元函数的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、如果绘出插值函数的曲线,观察其光滑性。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量与时间t的拟合曲线。t(分)051015202530354045505501.272.162.863.443.874.154.374.514.584.024.64二、要求1、用最小二乘法进行曲线拟合;2、近似解析表达式为;3、打印出拟合函数,并打印出与的误差,;4、另外选取一个近似表达式,尝试拟合效果的比较;5、*绘制出曲线拟合图。三、目的和意义1、掌握曲线拟合的最小二乘法;2、最小二乘法亦可用于解超定线代数方程组;3、探索拟合函数的选择与拟合精度间的关系。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验三数值积分与数值微分一、基本题选用复合梯形公式,复合Simpson公式,Romberg算法,计算(1)(2)(3)二、应用题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)用逐次分半辛普生法,并确定最佳步长。四、目的和意义1、深刻认识数值积分法的意义;2、明确数值积分精度与步长的关系;3、根据定积分的计算方法,结合专业考虑给出一个二重积分的计算问题。五、实验学时:2学时六、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验四线方程组的直接解法一、问题提出给出下列几个不同类型的线性方程组,请用适当算法计算其解。设线性方程组设对称正定阵系数阵线方程组三对角形线性方程组二、要求1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法;平方根法与改进平方根法;追赶法求解(选择其一);2、应用结构程序设计编出通用程序;3、比较计算结果,分析数值解误差的原因;4、尽可能利用相应模块输出系数矩阵的三角分解式。三、目的和意义1、通过该课题的实验,体会模块化结构程序设计方法的优点;2、运用所学的计算方法,解决各类线性方程组的直接算法;3、提高分析和解决问题的能力,做到学以致用;通过三对角形线性方程组的解法,体会稀疏线性方程组解法的特点。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验五解线性方程组的迭代法一、问题提出对实验四所列目的和意义的线性方程组,试分别选用Jacobi迭代法,Gauss-Seidel迭代法和SOR方法计算其解。二、要求1、体会迭代法求解线性方程组,并能与消去法做以比较;2、分别对不同精度要求,如由迭代次数体会该迭代法的收敛快慢;3、对方程组2,3使用SOR方法时,选取松弛因子ω=0.8,0.9,1,1.1,1.2等,试看对算法收敛性的影响,并能找出你所选用的松弛因子的最佳者;4、给出各种算法的设计程序和计算结果。三、目的和意义1、通过上机计算体会迭代法求解线性方程组的特点,并能和消去法比较;2、运用所学的迭代法算法,解决各类线性方程组,编出算法程序;3、体会上机计算时,终止步骤或k>(给予的迭代次数),对迭代法敛散性的意义;体会初始解,松弛因子的选取,对计算结果的影响。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验六非线性方程求根一、问题提出设方程有三个实根现采用下面六种不同计算格式,求f(x)=0的根或1、2、3、4、5、6、二、要求1、编制一个程序进行运算,最后打印出每种迭代格式的敛散情况;2、用事后误差估计来控制迭代次数,并且打印出迭代的次数;3、初始值的选取对迭代收敛有何影响;4、分析迭代收敛和发散的原因。三、目的和意义1、通过实验进一步了解方程求根的算法;2、认识选择计算格式的重要性;3、掌握迭代算法和精度控制;4、明确迭代收敛性与初值选取的关系。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验七矩阵特征值问题计算一、问题提出利用冪法或反冪法,求方阵的按模最大或按模最小特征值及其对应的特征向量。设矩阵A的特征分布为:且试求下列矩阵之一(1)求,及取结果(2)求及取结果:(3)求及取结果(4)取这是一个收敛很慢的例子,迭代次才达到结果(5)有一个近似特征值,试用幂法求对应的特征向量,并改进特征值(原点平移法)。取结果二、要求1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计;2、会用原点平移法改进算法,加速收敛;对矩阵B=A-PI取不同的P值,试求其效果;3、试取不同的初始向量,观察对结果的影响;4、对矩阵特征值的其它分布,如且如何计算。三、目的和意义1、求矩阵的部分特征值问题具有重要实际意义,如求矩阵谱半径,稳定性问题往往归于求矩阵按模最小特征值;2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧;3、问题中的题(5),反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征向量。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验八常微分方程初值问题数值解法一、基本题科学计算中经常遇到微分方程(组)初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题:(1)分别取h=0.1,0.2,0.4时数值解。初值问题的精确解。(2)用r=3的Adams显式和预-校式求解取步长h=0.1,用四阶标准R-K方法求值。(3)用改进Euler法或四阶标准R-K方法求解取步长0.01,计算数值解,参考结果。(4)利用四阶标准R-K方法求二阶方程初值问题的数值解(I)(II)(III)(IV)二、应用题1.小型火箭初始质量为900千克,其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力.当燃料用尽时引擎关闭。设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数为0.4(千克/米).重力加速度取9.8米/秒2.建立火箭升空过程的数学模型(微分方程);求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时间和高度.2.小型火箭初始质量为1200千克,其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生40000牛顿的恒定推力.当燃料用尽时引擎关闭。设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数记作k,火箭升空过程的数学模型为其中为火箭在时刻t的高度,m=1200-15t为火箭在时刻t的质量,T(=30000牛顿)为推力,g(=9.8米/秒2)为重力加速度,t1(=900/15=60秒)为引擎关闭时刻.今测得一组数据如下(t~时间(秒),x~高度(米),v~速度(米/秒)):t1011121314151617181920x10701270148017001910214023602600283030703310v190200210216225228231234239240246现有两种估计比例系数k的方法:1.用每一个数据(t,x,v)计算一个k的估计值(共11个),再用它们来估计k。2.用这组数据拟合一个k.请你分别用这两种方法给出k的估计值,对方法进行评价,并且回答,能否认为空气阻力系数k=0.5(说明理由).三、要求1、根据初值问题数值算法,分别选择二个初值问题编程计算;2、试分别取不同步长,考察某节点处数值解的误差变化情况;3、试用不同算法求解某初值问题,结果有何异常;4、分析各个算法的优缺点。四、目的和意义1、熟悉各种初值问题的算法,编出算法程序;2、明确各种算法的精度与所选步长有密切关系;3、通过计算更加了解各种算法的优越性。五、实验学时:2学时六、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.数值分析实验报告数值分析实验报告实验一函数插值方法报告一、问题提出对于给定的一元函数的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、如果绘出插值函数的曲线,观察其光滑性。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.求解如下:一、编写插值函数结构程序Lagrange插值多项式M文件:lagrange1.mfunction[A1,LN,L1,B1]=lagrange1(X,Y)m=length(X);LN=ones(m,m);fork=1:mx1=1;fori=1:mifk~=ix1=conv(x1,poly(X(i)))/(X(k)-X(i));endend L1(k,:)=x1;B1(k,:)=poly2sym(x1)endA1=Y*L1;LN=Y*B1(1)、求五次Lagrange多项式,在主显示区,输入五次Lagrange多项式程序:>>X=[0.40.550.650.800.951.05];>>Y=[0.410750.578150.696750.901.001.25382];>>[A1,LN,L1,B1]=lagrange1(X,Y)>>plot(X,A1);>>F=poly2sym(A1)运行后,输出五次Lagrange多项式的结果:A1=121.6264-422.7503572.5667-377.2549121.9718-15.0845F=(2139673480305281*x^5)/17592186044416-(1859275536318005*x^4)/4398046511104+(9836621836743*x^3)(414796119737013*x^2)/1099511627776+(2145751274873259*x)/17592186044416-1061478972867847/70368744177664拉格朗日插值多项式的图如下:(2)、试构造Lagrange多项式。结果为,在主显示区,输入程序:>>X=[1234567];>>Y=[0.3680.1350.0500.0180.0070.0020.001];>>[A1,LN,L1,B1]=lagrange1(X,Y)>>plot(X,A1);>>F=poly2sym(A1)运行后,输出结果的Lagrange多项式的结果:A1=0.0001-0.00160.0186-0.11750.4419-0.96830.9950F=(4304240283865561*x^6)/73786976294838206464-(7417128346304051*x^5)/4611686018427387904+(223*x^4)/12000-(2821*x^3)/24000+(994976512675275*x^2)/2251799813685248-(19367*x)/20000+199/200Lagrange多项式的图如下:二、计算函数值计算函数值的主程序:lagrangezhi.mfunction[y,R]=lagrangezhi(X,Y,x,M)n=length(X);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;q1=1.0;c1=1.0;forj=1:nifj~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;ends=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;(1)、计算、的值。在主显示区,输入程序:>>x=0.596;M=1;X=[0.4,0.55,0.65,0.80,0.95,1.05];>>Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382];>>[y,R]=lagrangezhi(X,Y,x,M)运行结果:y=0.6257R=2.2170e-008在主显示区,输入程序:>>x=0.99;M=1;X=[0.4,0.55,0.65,0.80,0.95,1.05];>>Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382];>>[y,R]=lagrangezhi(X,Y,x,M)运行结果:y=1.0542R=5.5901e-008(2)、计算、的值在主显示区,输入程序:>>x=1.8;M=1;X=[1,2,3,4,5,6,7];>>Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001];>>[y,R]=lagrangezhi(X,Y,x,M)运行结果:y=0.1648R=0.0059在主显示区,输入程序:>>x=6.15;M=1;X=[1,2,3,4,5,6,7];>>Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001];>>[y,R]=lagrangezhi(X,Y,x,M)运行结果:y=0.0013R=0.0042三、Newton插值多项式Newton插值多项式主程序M文件:Newton.mfunction[A,C,L,wcgs,Cw]=Newton(X,Y)n=length(X);A=zeros(n,n);A(:,1)=Y';s=0.0;p=1.0;q=1.0;c1=1.0;forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endb=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1;endC=A(n,n);b=poly(X(n));q1=conv(q1,b);fork=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C);Q=poly2sym(q1);symsMwcgs=M*Q/c1;Cw=q1/c1;在主显示区,输入的程序:>>x=[0.40.550.650.800.951.05];>>y=[0.410750.578150.696750.901.001.25382];>>[A,C,L,wcgs,Cw]=Newton(x,y)>>symsx;>>ezplot(L,[01.1]);运行结果如下,得到A=0.41080.57821.11600.69671.18600.28000.90001.35500.67600.99001.00000.6667-2.2944-7.4261-15.30201.25382.53827.486124.451463.7551121.6264C=121.6264-422.7503572.5667-377.2549121.9718-15.0845L=(8558693921221117*x^5)/70368744177664-(3718551072636019*x^4)/8796093022208+(5036350380412441*x^3)/8796093022208-(3318368957896111*x^2)/8796093022208+(536437818718315*x)/4398046511104-8491831782942691/562949953421312wcgs=(M*(x^6-(22*x^5)/5+(1583*x^4)/200-(3721*x^3)/500+(542206127247039*x^2)/140737488355328-(4682696525551953*x)/4503599627370496+4111390143022055/36028797018963968))/720Cw=0.0014-0.00610.0110-0.01030.0054-0.00140.0002牛顿插值多项式的图如下:在主显示区,输入的程序:>>x=[1234567];>>y=[0.3680.1350.0500.0180.0070.0020.001];>>[A,C,L,wcgs,Cw]=Newton(x,y)>>symsx;>>ezplot(L,[08]);运行结果如下,得到:A=0.36800.1350-0.23300.0500-0.08500.07400.0180-0.03200.0265-0.01580.0070-0.01100.0105-0.00530.00260.0020-0.00500.0030-0.00250.0007-0.00040.0010-0.00100.0020-0.00030.0005-0.00000.0001C=0.0001-0.00160.0186-0.11750.4419-0.96830.9950L=(8608480567731121*x^6)/147573952589676412928-(7417128346304047*x^5)/4611686018427387904+(223*x^4)/12000-(2821*x^3)/24000+(7959812101402191*x^2)/18014398509481984-(19367*x)/20000+199/200wcgs=(M*(x^7-28*x^6+322*x^5-1960*x^4+6769*x^3-13132*x^2+13068*x-5040))/5040Cw=0.0002-0.00560.0639-0.38891.3431-2.60562.5929-1.0000牛顿插值多项式的图如下四、结果讨论和分析通过上机实验,我对主要函数常见的插值方法有了更深的认识。在编辑拉格朗日插值、牛顿插值以及分段插值MATLAB程序过程中,深入了解了它们各自的算法,为我以后在工程中应用插值函数奠定了基础。根据得出的数据知道,与拉格朗日插值相比牛顿插值的计算量要小,程序设计方便,并且牛顿插值曲线相对要光滑些。分段三次埃尔米特插值与分段低次插值函数都有一致收敛性但是低次插值要求给出节点上的导数值,其光滑性较差,改进这种插值以克服其缺点要用三次样条插值。三次样条插值具有良好的收敛性和稳定性,又有二阶光滑度。实验二函数逼近与曲线拟合报告一、问题提出从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量与时间t的拟合曲线。t(分)051015202530354045505501.272.162.863.443.874.154.374.514.584.024.64二、要求1、用最小二乘法进行曲线拟合;2、近似解析表达式为;3、打印出拟合函数,并打印出与的误差,;4、另外选取一个近似表达式,尝试拟合效果的比较;5、*绘制出曲线拟合图。三、目的和意义1、掌握曲线拟合的最小二乘法;2、最小二乘法亦可用于解超定线代数方程组;3、探索拟合函数的选择与拟合精度间的关系。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.求解如下:实验步骤(一)算法流程构造a1、a2、a3的线性方程组------构造误差平方和------对a1、a2、a3求偏导数------令偏导为零求得a1、a2、a3的值。(二)编程步骤与分析1.绘制数据点(t,yi)的散点图输入程序为:t=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64].*1e-4plot(t,y,'r*'),legend('实验数据(t,yi)')xlabel('x'),ylabel('y'),title('数据点(t,yi)的散点图'),显示结果为:2.求参数a1、a2、a3的解析表达式计算在处的函数值,即输入程序symsa1a2a3t=[0510152025303540455055];fi=a1.*t+a2.*t.^2+a3.*t.^3运行后屏幕显示关于a1,a2,a3的线性方程组:fi=[0,5*a1+25*a2+125*a3,10*a1+100*a2+1000*a3,15*a1+225*a2+3375*a3,20*a1+400*a2+8000*a3,25*a1+625*a2+15625*a3,30*a1+900*a2+27000*a3,35*a1+1225*a2+42875*a3,40*a1+1600*a2+64000*a3,45*a1+2025*a2+91125*a3,50*a1+2500*a2+125000*a3,55*a1+3025*a2+166375*a3]构造误差平方和:y=[01.272.162.863.443.874.154.374.514.584.024.64].*1e-4;fy=fi-y;fy2=fy.^2;J=sum(fy.^2);运行后屏幕显示误差平方和如下:J=12650*a1^2+1089000*a1*a2+49967500*a1*a3-(7819978335372091569501*a1)/28823037615171174400000+24983750*a2^2+2386725000*a2*a3-(31331074233255294718193*a2)/2882303761517117440000+58593218750*a3^2-(274377591928296252150123*a3)/576460752303423488000+520374483464852566590953249225508026224249/332306998946228968225951765070086144000000000000为求a1,a2,a3使达到最小,只需利用极值的必要条件(k=1,2,3),得到关于a1,a2,a3的线性方程组,这可以由下面的MATLAB程序完成,即输入程序Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3),运行后屏幕显示J分别对a1,a2,a3的偏导数如下Ja11=25300*a1+1089000*a2+49967500*a3-27131/100000Ja21=1089000*a1+49967500*a2+2386725000*a3-217403/20000Ja31=49967500*a1+2386725000*a2+117186437500*a3-1903877/4000解线性方程组Ja11=0,Ja21=0,Ja31=0,Ja41=0,输入下列程序A=[25300,1089000,49967500;1089000,49967500,2386725000;49967500,2386725000,117186437500];B=[27131/100000,217403/20000,1903877/4000];C=B/A运行后求得系数C如下:C=2.6569e-005-5.2948e-0073.5168e-009故所求的拟合曲线为f(x)=2.6569×10-5t-5.2948×10-7t2+3.5168×10-9t33.估计其误差,并作出拟合曲线和数据的图形输入程序:t=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64].*1e-4;n=length(t);f=2.6569*1e-5*t-5.2948*1e-7*t.^2+3.5168*1e-9*t.^3;tx=0:1:100;F=2.6569*1e-5*tx-5.2948*1e-7*tx.^2+3.5168*1e-9*tx.^3;fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)plot(t,y,'r*'),holdon,plot(tx,F,'b-'),holdofflegend('数据点(t,y)','拟合曲线y=f(x)'),xlabel('x'),ylabel('y'),title('数据点(t,y)和拟合曲线y=f(x)的图形')运行后屏幕显示数据(t,y)与拟合函数f的最大误差Ew,平均误差E1和均方根误差E2及其数据点和拟合曲线y=f(x)的图形结果为:最大误差:Ew=4.2350e-005平均误差:E1=9.0785e-006均方根误差:E2=1.4684e-0054.用软件自带函数进行三次多项式拟合,比较拟合效果编写程序如下:t=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64].*1e-4;tx=0:1:100;F=2.6569*1e-5*tx-5.2948*1e-7*tx.^2+3.5168*1e-9*tx.^3;plot(t,y,'r*'),holdon,plot(tx,F,'b-')F1=polyfit(t,y,3)F1=F1(1)*tx.^3+F1(2)*tx.^2+F1(3)*tx+F1(4);plot(tx,F1,'+'),holdofflegend('数据点(t,y)','拟合曲线y=f(x)','拟合曲线y=f1(x)'),xlabel('x'),ylabel('y'),title('数据点(t,y)、拟合曲线y=f(x)、y=f1(x)图形')结果如下图:用软件拟合的三次多项式参数为F1(1)=3.4364e-009F1(2)=-5.2156e-007F1(3)=2.6340e-005F1(4)=1.7839e-006故软件拟合出的三次多项式函数为:F1(x)=2.6340×10-5t-5.2156×10-7t2+3.4364×10-9t3与求得给定表达式的拟合函数f(x)=2.6569×10-5t-5.2948×10-7t2+3.5168×10-9t3相比较,对应系数值非常接近,证明了所求解解析表达式的正确性。实验三数值积分与数值微分一、基本题选用复合梯形公式,复合Simpson公式,Romberg算法,计算(1)(2)(3)二、要求1、编制数值积分算法的程序;2、对基本题,分别取不同步长,试比较计算结果(如n=10,20等),并比较其结果;4、对应用题,用给定精度ε,试用(1)用逐次分半梯形法。(2)用逐次分半辛普生法,并确定最佳步长。三、目的和意义1、深刻认识数值积分法的意义;2、明确数值积分精度与步长的关系;3、根据定积分的计算方法,结合专业考虑给出一个二重积分的计算问题。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验过程:1、编写复合辛普森求积公式的M文件:functions=simpr1(f,a,b,n)%UNTITLEDSummaryofthisfunctiongoeshere%本程序为复合辛普生法求数值积分,f为被积函数,以内联形式定义该函数,即用inline()方式%a,b分别为积分的上下限%s为所求积分的结果%n为积分子区间数h=(b-a)/(2*n);s1=0;s2=0;fork=1:nx=a+h*(2*k-1);s1=s1+feval(f,x);endfork=1:(n-1)x=a+h*2*k;s2=s2+feval(f,x);ends=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;end在命令窗口输入以下程序:1)>>f=inline('exp(x)/(4+x^2)');>>s=simpr1(f,0,1,8)得出:s=1.58332)>>f='log(1+x)/(1+x^2)';>>s=simpr1('f',0,1,8)得出:s=1.58333)>>f='sqrt(4-(sin(x))^2)';>>s=simpr1('f',0,1/4,8)得出:s=0.2562编写龙贝格求积公式的M文件:function[R,quad,err,h]=romber(f,a,b,n,delta)%UNTITLED2Summaryofthisfunctiongoeshere%本程序为龙贝格求积程序,f为被积函数,其用内联函数形式定义%a,b分别为积分的上下限%delta为允许的误差,R为T数表,quad为求积结果M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(feval(f,a)+feval(f,b))/2;while((err>delta)&(J<n))|(J<4)J=J+1;h=h/2;s=0;forp=1:Mx=a+h*(2*p-1);s=s+feval(f,x);endR(J+1,1)=R(J,1)/2+h*s;M=2*M;fork=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1);enderr=abs(R(J,J)-R(J+1,k+1));endquad=R(J+1,J+1)end在命令窗口输入以下程序:1)>>f='exp(x)/(4+x^2)';>>[R,quad,err,h]=romber('f',0,1,8,1e-9)得出结果:quad=1.5833R=2.000000001.68751.58330001.60941.58331.5833001.58981.58331.58331.583301.58501.58331.58331.58331.5833quad=1.5833err=0h=0.0625>>f='sqrt(4-(sin(x))^2)';>>[R,quad,err,h]=romber('f',0,1/4,8,1e-9)得出quad=0.2562R=0.259800000.25710.25620000.25640.25620.2562000.25620.25620.25620.256200.25620.25620.25620.25620.2562quad=0.2562err=0h=0.0156>>f='log(1+x)/((1+x^2)';>>[R,quad,err,h]=romber('f',0,1,8,1e-9)得出:quad=1.5833R=2.000000001.68751.58330001.60941.58331.5833001.58981.58331.58331.583301.58501.58331.58331.58331.5833quad=1.5833err=0h=0.06253编写复合梯形求积公式的M文件:functionI=T(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);在命令窗口输入以下程序:1)>>x=0:0.1:1;>>y=(sin(x))./x;>>y(1)=1;>>I=T(x,y)运行可得:I=0.94582)>>x=0:0.1:1;>>y=exp(x)./(4+x.^2);>>I=T(x,y)运行可得:I=0.3909>>x=0:0.1:1;>>y=exp(x)./(4+x.^2);>>I=T(x,y)I=0.39083)>>x=0:0.1:1;>>y=log(1+x)./(1+x.^2);>>I=T(x,y)运行可得:I=0.2713>>x=0:0.01:1;>>y=log(1+x)./(1+x.^2);>>I=T(x,y)I=0.2722由以上结果可看出,若区间的份数不一样,则用同种方法算出的结果会有稍微的变化,但不会相差太大。>>x=0:0.01:1;>>y=sqrt(4-(sin(x)).^2);>>I=T(x,y)I=1.9298实验四线方程组的直接解法一、问题提出给出下列几个不同类型的线性方程组,请用适当算法计算其解。设线性方程组设对称正定阵系数阵线方程组三对角形线性方程组实验步骤:列主元高斯消去法的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_max<1e-10index=0;return;end%交换两行ifr>kforj=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;endendifmax1<1e-10falg='failure';return;end%交换两行ifr>kforj=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>(给予的迭代次数),对迭代法敛散性的意义;体会初始解,松弛因子的选取,对计算结果的影响。四、实验学时:2学时五、实验步骤:1.进入C或matlab开发环境;2.根据实验内容和要求编写程序;3.调试程序;4.运行程序;5.撰写报告,讨论分析实验结果.实验过程:编写雅克比迭代的M文件:functionx1=Jacobi(x0,a,b)e=1e-5;M=1000;n=length(x0);k=0;whilek<Mfori=1:ns=0;forj=1:nifj~=is=s+a(i,j)*x0(j);endendx1(i)=(b(i)-s)/a(i,i);endt=0;fori=1:nt=t+abs(x1(i)-x0(i));endift<ereturnendx0=x1;k=k-1;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];>>x0=[0.5-0.5011020.5-0.51];>>>>x=Jacobi(x0,a,b)运行可得:X=1.0000-1.00000.00001.00002.00000.00003.00001.0000-1.00002.0000实验六非线性方程求根一、问题提出设方程有三个实根现采用下面六种不同计算格式,求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)ifnargin<4itmax=100;endifnargin<3ep=1e-5;endindex=0;k=1;whilek<itmaxx1=x;x=feval(fun,x);ifabs(x-x1)<epindex=1;break;endk=k+1;endx_star=x;it=k;第二步:分别用上述六种形式的表达式计算方程的根,初值分别取2、0、-2结果如下。1、>>fun=inline('(3*x+1)^(1/3)');>>[x_star,index,it]=DD(fun,2)结果为:x_star=1.8794index=1it=92、>>fun=inline('(3*x+1)/x^2');>>[x_star,index,it]=DD(fun,-2)结果为:x_star=-1.5321index=1it=36>>[x_star,index,it]=DD(fun,0)结果为:x_star=NaNindex=0it=1002、>>fun=inline('(x+1)^(1/3)');>>[x_star,index,it]=DD(fun,-2)结果为:x_star=1.3247+0.0000iindex=1it=93、>>fun=inline('(x^3-1)/x');>>[x_star,index,it]=DD(fun,2)结果为:x_star=NaNindex=0it=1004、>>fun=inline('1/(x^2-3)');>>[x_star,index,it]=DD(fun,-1)结果为:x_star=-0.3473index=1it=65、>>fun=inline('(3+1/x)^1/2');>>[x_star,index,it]=DD(fun,-1)结果为:x_star=1.7808index=1it=96、>>fun=inline('x-(1/3)*((x^3-3*x-1)/(x^2-1))');>>[x_star,index,it]=DD(fun,2)结果为:x_star=1.8794index=1it=4>>[x_star,index,it]=DD(fun,0)结果为:x_star=-0.3473index=1it=4>>[x_star,index,it]=DD(fun,-2)结果为:x_star=-1.5321index=1it=5五、实验结论对于非线性方程,求它的解析解有时候是很困难的,但采用数值方法可以很容易地求它的近似解。此次实验就是采用迭代法求非线性方程的根。对于一个非线性方程,选用不同的迭代形式,因为其收敛程度不一样,造成其效率与精确度有很大的差别。实验七矩阵特征值问题计算一、问题提出利用冪法或反冪法,求方阵的按模最大或按模最小特征值及其对应的特征向量。设矩阵A的特征分布为:且试求下列矩阵之一(1)求,及取结果(2)求及取结果:(3)求及取结果(4)取这是一个收敛很慢的例子,迭代次才达到结果(5)有一个近似特征值,试用幂法求对应的特征向量,并改进特征值(原点平移法)。取结果二、要求1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计;2、会用原点平移法改进算法,加速收敛;对矩阵B=A-PI取不同的P值,试求其效果;3、试取不同的初始向量,观察对结果的影响;
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 2024年道路旅客运输从业资格证模拟考试题
- 2024年呼和浩特客运资格证都考些什么
- 专题05天气与气候-2024年中考地理专练(原卷版)
- 吉首大学《流体力学与液压传动》2021-2022学年第一学期期末试卷
- 吉首大学《电子商务概论》2021-2022学年第一学期期末试卷
- 《机加工艺方案设计与实施》考试卷A卷及答案
- 吉林艺术学院《影视声音基础》2021-2022学年第一学期期末试卷
- 吉林艺术学院《视觉特效制作与合成》2021-2022学年第一学期期末试卷
- 转让个人板车协议书范本模板
- 村民占地调节协议书范文范本
- 北京市西城区2023-2024学年高一下学期期末英语试题(解析版)
- 金华市金投集团有限公司招聘笔试题库2024
- 中国中煤笔试
- 人教版pep五上《Unit 4 What can you do》说课稿
- 4.2+在实践中追求和发展真理 课件 高中政治统编版 必修四 哲学与文化
- Chat GPT 科普知识讲解
- 山西退役军人事务厅事业单位笔试真题2024
- DBJ50-T-271-2017 城市轨道交通结构检测监测技术标准
- 医学美容技术专业《医学美学导论》课程标准
- 第二单元缤纷舞曲 主题教学设计 2023-2024学年人音版七年级上册教案1000字
- 幼儿园立德树人实施方案及措施(2篇)
评论
0/150
提交评论