还剩23页未读，继续阅读
本资源只提供10页预览，全部文档请下载后查看！喜欢就下载吧，查找使用更方便
文本内容:
《数值分析》实验报告姓名______________学号专业指导教师刘建生教授日期年月日20151225实验_插值Lagrange/newton=于一:对于给定的一元函数的个节点值，几试用公式求y=/x n+1yj x0j=0,1,Lagrange其插值多项式或分段二次插值多项式Lagrange数据如下求五次
0.
40.
550.
650.
800.
951.
050.
410750.
578150.
696750.
901.
001.25382计算的值提示结果为了/
0.596,/
0.
990.
5960.625732,/
0.99«
1.
0542312345670.
3680.
1350.
0500.
0180.
0070.
0020.001试构造多项式，计算的了值提示结果为Lagrange Le
1.8,/
6.15/L8a
0.164762,f
6.15«
0.001266二实验程序及注释程序MATLAB functionf=lagrange x0,yO,xn=lengthxO;m=lengthyO;format longs=
0.0;for k=l:nP=
1.0;for j=l:n ifj~=ks=0;for i=l:nx=a+h*2*i-l;s=s+zhuci_tixing_fx;endT=T0/2+h*s;n=2*n;delt=T-T0;endjifen=T;num=n;end执行程序后的结果实验分析逐次分半的积分算法具有很高的精度，对于复杂的工程实践问题具有很高的应用价值同时，在利用复合梯形公式、复合公式以及算法等计算数值积分时，精度已经Simpson Romberg很高，并且随着步长越小，积分精度越高另外，当给定一个精度的时候，我们还可以利用自适应步长法，确定所需要的最佳步长和结果同样在对比中可见公式的收敛速度比梯形公式的收敛速度快simpson实验四线性方程组的直接解法
一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解、设线性方程组1参考解:X*=1,—1,0,1,2,0,3,1,—1,
27、设对称正定阵系数阵线方程组2参考解:£=—1,0,21,—1,0,2,三对角形线性方程组3,参考解:，2,1,-3,0,1,-2,3,0,1,-1要求、对上述三个方程组分别利用顺序消去法与列主元消去法；平方根法与1Gauss Gauss改进平方根法；追赶法求解选择其一；、应用结构程序设计编出通用程序；
2、比较计算结果，分析数值解误差的原因；
3、尽可能利用相应模块输出系数矩阵的三角分解式4
二、问题分析直接法就是经过有限步算术运算，可求得线性方程组精确解得方法若计算过程中没有舍入误差但是实际运算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解我们经常采用的一些方法，如顺序消去法与列主元消去法、平方根法与改进Gauss Gauss平方根法、追赶法等等
三、实验程序及注释%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%列主兀消去法Gaussfunction[x,det,index]=GaussA,b%求线性方程的列主元消去法Gauss为方程组的系数矩阵%A为方程组的右端项%b为方程组的解%x为系数矩阵的行列式的值%det A为指标变量，表示计算失败，表示计算成功%index index=0index=l为系数矩阵的行列式的值%det A[n,m]=sizeA;nb=lengthb;%当方程组行与列的维数不相等时，停止计算，并输出错误的信息if n〜=m的行与列必须相等！，eiroifAreturn;end%当方程组与右端项的维数不匹配时，停止计算，输出错误信息if m〜=nbeirorthe columnsof Amust beequal thelength ofb;return;end%开始计算，先赋初值index=1;det=1;x=zerosn,l;for k=l:n-l%选主元a_max=0;for i=k:nif absAi,ka_maxa_max=absAi,k;r=i;endendif a_max le-10index=O;return;end%交换两行if rk;for j=k:nz=Ak,j;Ak,j=Arj;Ar,j=z;endz=bk;bk=br;br=z;det=-det;end%消元过程for i=k+l:nm=Ai,k/Ak,k;for j=k+l:nAi,j=Ai,j-m*Ak,j;endbi=bi-m*bk;enddet=det*Ak,k;enddet=det*An,n;%回代过程if absAn,nle-10index=O;return;endfor k=n:-l:lfor j=k+l:nbk=bk-Ak,j*xj；endxk=bk/Ak,k;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%改进平方根法function[L,D,index]=LDL_DecomA%求正定对称矩阵的改进平方根分解，也称分解A LDLT为要分解的矩阵%A为下三角阵%L%口为对角阵为指标变量，表示计算失败，表示计算成功%index index=0index=l n=length A;L=eyen;D=zerosn;d=zeros1,n;T=zerosn;index=1;for k=l:ndk=Ak,k;for j=l:k-ldk=dk-Lk,j*Tk,j;endif absdkle-10index=O;return;endfor i=k+l:nTi,k=Ai,k;for j=l:k-lTi,k=Ti,k-Ti,j*LkJ;endLi,k=Ti,k/dk;endendD=diagd;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%追赶法function x=followupA,b%采用追赶法求解%系数矩阵A%常数向量b%线性方程组的解xn=rankA;for i=l:nifAi,i==0对角元素为disperror:O1;return;endend;d=onesn/；a=onesn-l,l;c=onesn-l;fori=l:n-lai1=Ai+l,i;ci,l=Ai,i+l;di,l=Ai,i;enddn/=An,n;fori=2:ndi1=di,1-ai-11/di-1』*ci-1』；bi,1bi,1-ai-1,1/di-1,1^bi-1,1;endxn,l=bn,l/dn,l;fori=n-l:-l:l』*』/xi,l=bi/-ci xi+l di,l;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%第一问A-[42-3-1210000;86-5-3650100;42-2-132-103l;0-215-13-1194;-426-167-3323;86-8571726-35;02-13-425301;1610-11-917342-122;462-713920124;00-l8-3-24-863-1]b=[51232346133819-21][x,det,index]=GaussA,b[L,U,P]=luA第二问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]1[L,D,index]=LDL_DecomAR=L*bm=D*Lm\r[L,U,P]=luA第三问A=[4-100000000;-l4-10000000;0-14-1000000;00-14-100000;000-14-10000;0000-14-1000;00000-14-100;000000-14-10;0000000-14-l;00000000-14]b=[75-1326-1214-45-5],x=followupA,b[L,U,P]=luA
四、实验数据结果及分析、第一问中的方程组属于一般的线性方程组，这里采用列主元高斯消去法，求得对比精1确值』,-、我们发现差距很大估计是所给的参考值错误同样利用函数x*=l,-1,0,21,0,2ki可以系数矩阵的三角分解式A、该方程组为三对角形线性方程组，可以采用追赶法解决类似的问题得到的结果2对比题中的参考值」,—丁，精度很高三角分解为x*=2,1,—3,0,1,—2,3,01
五、实验结论本实验运用最基本的算法高斯消去法以及某些变形这类方法是解低阶稠密矩阵方程组以及某些大型稀疏矩阵方程组的有效方法列主元消去法比顺序消去法更加能保证结果精度，因为列主元消去法始终保持除Gauss数最大，避免了因为计算机内存结构而引起的“大数”吃“小数”，故在工程中比较常用实验五解线性方程组的迭代法一问题提出对实验四所列目的和意义的线性方程组，试分别选用迭代法，Jacobi Gauss-Seidel迭代法和方法计算其解SOR实验程序代码.迭代法计算线性方程组%jacobifunction[x k]=j acobi_fA,b,xO,tol maxl=300;rD=diagdiagA;L=-trilA,-1;U=-triuA,1;B=D\L+U;f=D\b;x=B*xO+f;k=l;while normx-xO=tolxO=x;x=B*xO+f;k=k4-l;ifk=maxl•迭代次数超过次，方程组可能无解，；disp300return;end二，％高斯一一赛德尔迭代法计算线性方程组function[x k]=guass_seid_fA,b,xO,tolzmaxl=300;D=diagdiagA;L=-trilA,-1;U=-triuA,1;G=D-L\U;f=D-L\b;x=G*xO+f;k=l;while normx-xO=tolx0=x;x=G*xO+f;k=k+l;ifk=maxl迭代次数太多，可能不收敛disp Dend%[k,xf]；end三，％超松弛迭代法方法函数sor%超松弛迭代法的统一方程function[x,k]=sor_fA,b,xO,w,tolmax=300;if w=0||w=2error;return;endD=diagdiagA;L=-trilA,-1;U=-triuA,1;B=invD-L*w*1-w*D+w*U;f=w*inv D-w*L*b;x=B*xO+f;k=l;while normx-xO=tolxO=x;x=B*xO+f;k=k+l;if kmax「迭代次数过多，方法可能不收敛disp sorD;return;end%[k,xf];%输出每一步迭代的结果，可不要End实验结果对比题中的参考值精度很高X*=2,1,-3,0,1,-2,3,0,1,—1\实验六非线性方程求根一问题提出设方程有三个实根；芯=/x=/—3x_1=0X=
1.8793,—
0.34727,芯现采用下面六种不同计算格式，求的根或芯=—
1.53209fx=O X3x+l1X=IX实验程序代码本程序采用牛顿法求解非线性方程组也可自行编写简单迭代算法进行求解，此处没有给出程序代码如下迭代法求解非线性方程%newton function[gen,time]=newtonf,xO,tol if nargin==2tol=l,Oe-5;enddf=difff;xl=xO;time=0;wucha=
0.1;whilewuchatol time=time+l;fx=subsf,xl;df=subsdf,xl;gen=xl-fx/df;wucha=absgen-xl;xl=gen;endend分别带入各个迭代格式得由收敛准则知道，当且仅当迭代方程满足以上的收敛条件时，非线性方程组才有^xLl解实验七矩阵特征值问题计算
一、问题提出利用幕法或反幕法，求方阵旬…的按模最大或按模最小特征值及其对应的4=特征向量设矩阵的特征分布为A国一氏以阖之…之％同且例=3/M试求下列矩阵之一---121求，及西1A=2—41411-6结果为
46.42106,“—
0.046152,—
0.374908,14-273—81-25117417252A=31261-15543231754取⑼，一=1,0,1,0,0,1£=105结果~2-1-12—13A=-1求及王2-14-12-1-12取44=LLUMgcr结果；I
3.73214取J=U,1,17£=10-2丁结果p1,
2.501460-
0.757730-
2.5642124«-
8.028578355这是一个收敛很慢的例子，迭代次才达到120021-411—6」有一个近似特征值—
6.42,试用基法求对应的特征向量,并改进特征值原点平移法取J°=1,UT£=17结果；乏I-
6.42107-
0.0461465-
0.37918运用指令MATLAB；也可用幕法求解,[V,D]=eigA程序在下面给出，可自行调试进行求解:123下面给出幕法的算法程序已调试function[namda,q,time,date_na,date_q]=mifaA,tol,xO if nargin==1tol=le-7;temp=lengthA;xO=onestemp,1;endif nargin==2temp=lengthA;xO=onestemp,1;endml=0;u=xO;time=0;while time500v=A*u;[vmax,i]=maxabsv;m=v i;u=v/m;if absm-mltol break;endml=m;time=time+l;date_natime,1=m;date_q:,time=u;endnamda=m;q=u;end实验八常微分方程初值问题数值解法问题的提出1科学计算中经常遇到微分方程（组）初值问题，需耍利用法，改进法（预Euler Euler Euler测一校正法），四阶方法求其数值解，诸如以下问题Rung-Kutta（）；y0=30x2分别取时数值解初值问题的精确解h=y=j4+5e3p=p*x-xOj/xOk-xO j;endends=s+yOk*p;Endf=s;end结果运行:»xO=[
0.
40.
550.
650.
800.
951.05];yO=[
0.
410750.
578150.
696750.
901.
001.25382];x=
0.596:lagrange xO^yO,xans=
0.625732377952666»xO=[
0.
40.
550.
650.
800.
951.05];yO=[
0.
410750.
578150.
696750.
901.
001.25382]x=
0.99:lagrange xO,yO,x ans=
1.054229770812718结果与提示值完全吻合，说明插值多项式的精度是很Lagrange的;WI J同时，若采用三点插值和两点插值的方法，用三点插值的精度更高若同时采用两点插值，选取的节点距离越近，精度越高x三采用插值进行计算newton算法程序如下format long;x0=[
0.
40.
550.
650.
800.
951.05];yO=[O.
410750.
578150.
696750.
901.
001.25382];x=
0.596;n=maxsizexO;y=yO1;%dispy;s=l;dx=yO;for i=l:n-
1、问题的分析
2、利用公式1Euler当取定初值时，即可由上式递推计算出y0=3yl,y2,...,yn.、利用改进的公式2Euler利用四阶公式3Rung-Kutta、实验主程序及部分注释3法编程主程序
3.1Euler法%%Eul㊀rfunction[x,y]=Eulerfun,span,yO,h x=span1:h:span2;y1=y0;for n=l:length x-1yn+1=yn+h*fevalfun,xn,yn;endR=[x Ty*]改进法编程主程序
3.2Euler%%改进的法预测一校正法EulerEulerfunction[x,y]=Euler1fun,span,yO h x=spanf1:h:span2;y1=y0;for n=l:lengthx-1kl=fevalfun,xn,yn;女通吗好理砒短璟*+h*kl;y n+1=yn+h*kl+k2/2;h_㊀nd y〃+]=%+-[/%„,M+“£+]，算+]]%%x=x，;y=yf「，一・・・D”0,1,2,R=[x Ty]end四阶法编程主程序
3.3Rung-Kutta%%标准四阶龙格库塔法function[x y]=RKfun,xspan,yO,hx=xspan1:h:xspan2y1=y0;zfor n=l:lengthx-1kl=fevalfun,xn,yn;k2=fevalfun,xn+h/2,yn+h/2*kl;k3=fevalfun,xn+h/2,y n+h/2*k2;k4=fevalfun,xn+h,yn+h*k3;yn+l=yn+h*kl+2*k2+2*k3+k4/6;endx=x,;y=y1;绘制出积分曲线图
3.4图-法与精确值比较1Euler图改进法与精确值比较-2Euler图法与精确值比较-3R-K运行程序
3.5法运行程序
3.
5.1Euler在的中键入如下内容:Matlab CommandWindow»format long»format long»format longfun=inline4*x/y-x*y;fun=inline4*x/y-x*y;fun=inline4*x/y-x*y;[x,y]=Eulerfun,[0,2],3,
0.1[x,y]=Eulerfun,[x,y]=Eulerfun,[0,2],3,
0.4回车结果xl yh=
0.1[0,x2]2,3,
0.2yh=
0.2x3yh=
0.4改进法运行程序
3.
5.2Euler在的中键入如下内容:Matlab CommandWindow»format long»format long»format longfun=inline4*x/y-x*y;fun=inline4*x/y-x*y;fun=inline4*x/y-x*y;[x,y]=Eulerlfun,[0,2],3,
0.1[x,y]=Euler1fun,[0,2],3,
0.2[x,y]=Euler1fun,[0,2],3,
0.4回车结果xl yh=
0.1x2yh=
0.2x3yh=
0.4四阶法运行程序
3.
5.3Rung-Kutta在的中键入如下内容Matlab CommandWindow»format long»format long»format longfun=inline4*x/y-x*y;二fun=inline4*x/y-x*y;fun inline4*x/y-x*y;[x,y]=RKfun,[0,2],3,
0.4[x,y]=RKfun,[0,2],3,
0.1[x,y]=RKfun,[0,2],3,
0.2回车结果xl yh=
0.1x2yh=
0.2x3yh=
0.4dxO=dx;for j=l:n-idxj=dxOj+1-dxOj/xOi+j-xOj;enddf=dx1;s=s*x-xO i;%计算y=y+s*df;%%dispy;enddispy运行结果绘制出曲线图与结果相吻合所以法和法的思想是一样的适合理论分析，但newton LagrangeLagrange法不如法灵活如果节点个数改变，算法需要重Lagrange newtonLagrange新编写，而法克服这一缺点，所以应用更为灵活Newton实验二函数逼近与曲线拟合
一、问题提出在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间的拟合曲线t分钟K0510152025303540455055丫°
01.
272.
162.
863.
443.
874.
154.
374.
514.
584.
024.
64、用最小二乘法进行曲线拟合；
12、近似解析表达式为fx=at+a12+a3t3；
2、计算出拟合函数并列出出与的误差；3fx,fx yx、另外选取一个近似表达式，尝试拟合效果的比较；
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位精度小数format long%14plot（t,y；*r）%绘制被拟合数据点的离散图hold on（）绘制三次拟合函数图（其中是拟合之后的数据）plot t,yl,b%yl（时间（分钟））%注释轴xlabel tx含碳量八）注释轴ylabelC/10-4%y三次拟合图％注释图名title,9坐标系网格化grid%四次拟合程序（最小二乘法）（）调用最小二乘法的程序进行四次拟合并给出误差分析[p,s]=polyfit t,y,4%MATLAB位精度小数format long%14plot（t,y*『）％绘制被拟合数据点的离散图hold on,）绘制三次拟合函数图（其中是拟合之后的数据）plot y2,b%y2时间（分钟））％注释轴xlabelf tx含碳量八）注释轴ylabelC/104%y出贝四次拟合图）％注释图名1坐标系网格化grid%
四、实验数据结果及分析三次拟合可以得到其拟合多项式为=0・000036t3-
0.001556t+
0.739853t+
0.083883拟合函数与被拟合函数图之间的对比如下（）红色星号为原始数据；1（）带圈的曲线为最小二乘后而成的结果曲线2由此可见拟合函数与原函数离散数据点拟合成程度相当好，通过（）对拟合误[p,s]=polyfit t,y,n差进行分析，如图图2-2由此可知，三次拟合精度较好为了提高结果的可信度，我们另外选取一个近似表达式，尝试拟合效果的比较于是,进行四次拟合其中，拟合得到的多项式为%=
0.000000t4-
0.000092t3-
0.00697712+
0.4944321+
0.094872拟合如图
2.3图2-3同样对四次拟合进行误差分析可得图2-4由此可见，四次拟合误差（三次拟合误差），精度更高
0.
494930.50824
五、实验结论在用高阶多项式对某一函数进行曲线拟合时,•并不是拟合出来的多项式与被拟合函数在整个区间上都能符合,只能保证在输入数据所能达到的区间上及其附近.求得的多项式polyfit x可以最大限度在逼近原函数利用最小二乘法对本问题进行的曲线拟合精度较高，而且，在一般情况下，拟合的多项式次数越多，精度越高实验三数值积分与数值微分
一、问题提出计算下列积分值21I=j V4-sin^/«
1.5343916o1・八21=0=
10.9460831i%----------3I=f-dx_r lnl+x,z x，41=―—dxi要求、编制数值积分算法的程序；
1、分别用两种算法计算同一个积分，并比较其结果；
2、分别取不同步长=/〃，试比较计算结果如等；3S—n=10,
20、给定精度要求£，试用变步长算法，确定最佳步长
二、问题分析4由上可知这四个积分找不到用初等函数表示的原函数，直接计算起来很困难，因此我们考虑利用函数在若干点得函数值，近似地计算该函数在一个区间上得定积分这里采用的方法有三种复合梯形公式，复合公式，算法Simpson Romberg
三、实验程序及注释、复合梯形公式程序1MATLABfunction I=T_quadx,y%复化梯形求积公式，其中，%*为向量，被积函数自变量的等距节点；%丫为向量，被积函数节点处的函数值；n=lengthx;m=lengthy;ifn〜=merrorthe lengthof Xand Ymust beequal!*;return;endh=xn-xl/n-l;a=[l2*onesl,n-21];I=h/2*suma.*y;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%、复合公式程序2Simpson MATLABfunction I=S_quadx,y%*为向量，被积函数自变量的等距节点；%丫为向量，被积函数节点处的函数值；n=lengthx;m=lengthy;ifn〜=merrorCthe lengthof Xand Ymust beequal!*;return;end如果不能被整除，则调用复化梯形公式if remn-l,2〜=0%n-12I=T_quadx,y;return;endN=n-1/2;h=xn-x1/N;a=zerosl,n;for k=l:Na2*k-l=a2*k・l+l;a2*k=a2*k+4;a2*k+l=a2*k+l+l;end二I h/6*suma.*y;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%、算法程序3Romberg MATLABfunctionI=R_quad_iterfun,a,b,ep求积公式，其中，%Romberg为被积函数；%fun为积分区间端点，要求%a,b ab;%叩精度系数，缺省值为le-
5.if nargin4ep=le-5;endm=l;h=b-a;I=h/2*fevalfun,a+fevalfun,b;Tl,1=1;while1八N=2m-l;h=h/2;I=I/2;for i=l:N;I=I+h*fevalfun,a+2*i-1*h;endTm+l1=I;M=2*N;k=l;while M1;Tm+1,k+1=4八k*Tm+1,k-Tm,k/4Ak-l;M=M/2;k=k+l;endif absTk,k-Tk-l,k-lepbreak;endm=m+1;endI=Tk,k;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%、自适应步长梯形公式4functionI=R_quad_iterfun,a,b,ep%梯形递推求积公式，其中，为被积函数；%fun为积分区间端点，要求%a,b avb;精度系数，缺省值为%ep le-
5.if nargin4ep=le-5;end;N=l;h=b-a;T=h/2*fevalfun,a+fevalfun,b;while1h=h/2;I=T/2;for k=l:N;I=I+h*fevalfun,a+2*k-1*h;endif absI-Tep break;endN=2*N;T=I;end对四个积分求解1/4__________%对1=，j4-sin2%上/PI.5343916求解ox=0:1/40:1/4;y=sqrt4-sinx.A2;format long%调用复化梯形公式求得积分值I=T_quadx,y I=
0.652029调用复化辛普森公式求得积分值I=S_quadx,y%I=
0.844568fun=inlinesqrt4-sinx.A2,;调用龙贝格算法可得积分值I=R_quad_iterfun,0,1%1=
0.4987x=0:1/400:1/4;y=sqrt4-sinx.A2;%调用复化梯形公式求得积分值I=T_quadx,y I=
0.466684调用复化辛普森公式求得积分值I=S_quadx,y%I=
0.757532%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1・%对求解x/o=1,/p
0.9460831o xx=
0.1:
0.1:l;%可以得到值，由题可知y=sinx./x yf0=loy=[
10.
6468280.
3975310.
5537800.
5771630.
7208410.
8991730.
1768130.
3624400.
2919430.480790];x=0:
0.1:l;调用复化梯形公式可以求得I=T_quadx,y%1=
0.186691调用复化辛普森公式求得积分值I=S_quadx,y%1=
0.883807使用不同的步长x=
0.01:
0.01:l;%xy=sinx./x调用复化梯形公式可以求得I=T_quadx,y%1=
0.887360%调用复化辛普森公式求得积分值I=S_quadx,y1=
0.833507%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1%对浊求解1=1——x=0:
0.1:l;y=expx./4+x.A2;调用复化梯形公式可求得积分值I=T_quadx,y%I=
0.274504调用复化辛普森公式求得积分值I=S_quadx,y%1=
0.686445使用不同的步长x=0:
0.01:l;%xy=expx./4+x.A2;调用复化梯形公式可求得积分值I=T_quadx,y%1=
0.557538调用复化辛普森公式求得积分值I=S_quadx,y%1=
0.557538%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对汁之求积分:\1+/x=0:
0.1:l;y=logl+x./l+x.A2;调用复化梯形公式可求得积分值I=T_quadx,y%I=
0.750865调用复化辛普森公式求得积分值I=S_quadx,y%1=
0.573417使用不同的步长x=0:
0.01:l;%xy=logl+x./l+x.A2;调用复化梯形公式可求得积分值I=T_quadx,y%1=
0.272调用复化辛普森公式求得积分值I=S_quadx,y%I=
0.272%给定精度要求心试用变步长算法，确定最佳步长fun=inlinesqrt4-sinx.A2;I=R_quad_iterfun,0,1/4
二、应用题文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，
1.在两坐标轴上取天文测量单位一天文单位为地球到太阳的平均距离万里在五个不同的时间9300对小行星作了五次观察，测得轨道上五个点的坐标数据如下表所示P P PPP12345坐标X
5.
7646.
2866.
7597.
1687.408坐标y
0.
6481.
2021.
8232.
5267.408由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为现需要建立椭圆的方程以供研究分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程1组AX=bo用求低价方程组的指令求出待定系数%•电巧，2MATLAB A\b4%卫星轨道是一个椭圆，其周长的计算公式是3式中，是椭圆的半长轴，尺.”+方/石是地球中心与轨道中心椭圆中心的a°-22距离，一方其中为近地点距离，为远地点距离，为地球半径C=“”2h HR=6371km有一颗人造卫星近地点距离远地点距离试分别按下列方案计算卫星轨h=439km,H=2384km道的周长，误差限取为10
三、要求、编制数值积分算法的程序；
1、对基本题，分别取不同步长〃=〃—〃/〃，试比较计算结果如2n=10,20等，并比较其结果；、对应用题，用给定精度£，试用⑴用逐次分半梯形法用逐次分半辛普生法，并确定最佳42步长逐次梯形法的主函数程序function[j ifen,num]=zhuci_tixinga,b,tolif nargin==3eps=l.Oe-3;enddelt=l;n=l;h=b-a;T=h*zhuci_tixing_fa+zhuci_tixing_fb/2;while delt3*tolh=h/2;T0=T;。
个人认证
优秀文档
获得点赞 0
