2021年数值计算B大作业

2021年数值计算B大作业2021年数值计算B大作业PAGE/NUMPAGES2021年数值计算B大作业课程 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 课程名称:数值计算B设计题目:数值计算B大作业学号:姓名:完成时间:题目一:多项式插值某气象观察站在8:00（AM）开始每隔10分钟对天气作以下观察,用三次多项式插值函数（Newton）迫近以下曲线,插值节点数据如上表,并求出9点30分该地域温度（x=10）。x12345678y22.523.324.421.7025.228.524.825.4二、数学原理假设有n+1个不一样节点及函数在节点上值（x,y）,……（x,y）,插值多项式有以下形式:（1）其中系数（i=0,1,2……n）为特定系数,可由插值样条（i=0,1,2……n）确定。依据均差定义,把x看成[a,b]上一点,可得f(x)=f（）+f[]（）f[x,]=f[]+f[x,]（）……f[x,,…x]=f[x,,…x]+f[x,,…x](x-x)综合以上式子,把后一式代入前一式,可得到:f(x)=f[]+f[]（）+f[]（）（）+…+f[x,,…x]（）…(x-x)+f[x,,…x,]=N（x）+其中N（x）=f[]+f[]（）+f[]（）（）+…+f[x,,…x]（）…(x-x)（2）=f(x)-N（x）=f[x,,…x,]（3）=（）…(x-x)Newton插值系数（i=0,1,2……n）能够用差商表示。通常有[]（k=0,1,2,……,n）（4）把（4）代入（1）得到满足插值条件N（i=0,1,2,……n）n次Newton插值多项式N（x）=f（）+f[]（）+f[]（）（）+……+f[]（）（）…（）.其中插值余项为:介于之间。三、程序设计function[y,A,C,L]=newdscg(X,Y,x,M)%y为对应x值,A为差商表,C为多项式系数,L为多项式%X为给定节点,Y为节点值,x为待求节点n=length(X);m=length(x);%n为X长度fort=1:mz=x(t);A=zeros(n,n);A(:,1)=Y';s=0.0;p=1.0;q1=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));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));fork=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)=polyval(C,z);%输出y值endL(k,:)=poly2sym(C);%输出多项式>>symsM,X=[1,3,5,7];Y=[22.5,24.4,25.2,24.8];x=10;>>[y,A,C,L]=newdscg(X,Y,x,M)y=21.7313A=22.500000024.40000.95000025.0.4000-0.1375024.8000-0.-0.1500-0.0021C=-0.0021-0.11871.452121.1688L=-x^3/480-(19*x^2)/160+(697*x)/480+3387/160结果分析和讨论对于不超出三次插值多项式,x假如选择1,3,5,7这三个点能够得到很好三次插值多项式L=-0.0021x^3-0.1187x^2+1.4521x+21.1688。当x=10时,也即9点30分时温度为21.7317度,结果分析知此值应是偏小。对于选择不一样插值节点,能够得到不一样插值多项式,误差也不尽相同。完成题目体会与收获牛顿插值法关键一点就是对插值节点选择,经过本题编程很好加深了对其概念了解以及提升了应用牛顿插值法能力,学会了利用Matlab软件对牛顿插值法相关问题进行编程求解,对Matlab计算方法与程序编辑愈加熟悉。使我对这类问题了解有了很大提升。题目二:曲线拟合在某钢铁厂冶炼过程中,依据统计数据含碳量与时间关系,试用最小二乘法拟合含碳量与时间t拟合曲线,并绘制曲线拟合图形。t(分)051015202530354045505501.272.162.863.443.874.154.374.514.584.024.64数学原理从整体上考虑近似函数同所给数据点(i=0,1,…,m)误差(i=0,1,…,m)大小,常见方法有以下三种:一是误差(i=0,1,…,m)绝对值最大值,即误差向量∞—范数;二是误差绝对值和,即误差向量r1—范数;三是误差平方和算术平方根,即误差向量r2—范数;前两种方法简单、自然,但不便于微分运算,后一个方法相当于考虑2—范数平方,所以在曲线拟合中常采取误差平方和来度量误差(i=0,1,…,m)整体大小。数据拟合具体作法是:对给定数据(i=0,1,…,m),在取定函数类中,求,使误差(i=0,1,…,m)平方和最小,即=从几何意义上讲,就是寻求与给定点(i=0,1,…,m)距离平方和为最小曲线。函数称为拟合函数或最小二乘解,求拟合函数方法称为曲线拟合最小二乘法。在曲线拟合中,函数类可有不一样选择方法。三、程序设计>>x=[0,5,10,15,20,25,30,35,40,45,50,55];>>y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64]*10^(-4);>>p=polyfit(x,y,2)p=1.0e-04*-0.00240.20370.2305>>plot(x,y,'r')结果分析和讨论经过最小二乘法得到曲线拟合多项式是p=(-0.0024x^2+0.2037x+0.2305)*10-4。利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数,因为数据较少得到拟合曲线不够光滑,可能会存在一定误差。拟合曲线总体展现增函数趋势,与数据较为吻合。完成题目体会与收获曲线拟合较常见方法之一就包含最小二乘法,所以能够熟练使用Matlab进行数据曲线拟合变得尤为关键。经过完成本题拟合过程,对于最小二乘法曲线拟合操作愈加熟练,能够很好完成各类数据拟合。题目三:线性方程组求解分别利用LU分解;平方根法与改善平方根法;追赶法求解下列多个不一样类型线性方程组,并与正确值比较结果,分析误差产生原因。（1）设线性方程组二、数学原理将A分解为一个下三角矩阵L和一个上三角矩阵U,即:A=LU,其中L=,U=解Ax=b问题就等价于要求解两个三角形方程组:　⑴Ly=b,求y;  ⑵Ux=y,求x。设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵。L,U元素能够有n步直接计算定出。用直接三角分解法解Ax=b（要求A全部次序主子式都不为零）计算公式:①,,i=2,3,…,n.计算U第r行,L第r列元素(i=2,3,…,n): ②  ,i=r,r+1,…,n;③,i=r+1,…,n,且r≠n.求解Ly=b,Ux=y计算公式;④⑤三、程序设计functionf1=LUresolve(a,b);[n,n]=size(a);%行数z=size(b)%b行数l=[];%l矩阵u=[];%u矩阵forj=1:1:nu(1,j)=a(1,j);endfori=2:1:nl(i,1)=a(i,1)/u(1,1);endfori=2:1:(n-1)sum1=0;fork=1:1:(i-1)sum1=l(i,k)*u(k,i)+sum1;endu(i,i)=a(i,i)-sum1;sum2=0;sum3=0;forj=(i+1):1:nfork=1:1:(i-1)sum2=sum2+l(i,k)*u(k,j);sum3=sum3+l(j,k)*u(k,i);endu(i,j)=a(i,j)-sum2;l(j,i)=(a(j,i)-sum3)/u(i,i);endendfori=1:1:(n-1)l(i,i)=1;l(i,n)=0;endl(n,n)=1;sum4=0;fork=1:1:(n-1)sum4=sum4+l(n,k)*u(k,n);endu(n,n)=a(n,n)-sum4;l%输出l矩阵u%输出u矩阵y=l\b%输出向量yx=u\y%输出向量x>>a=[42-3-121000086-5-365010042-2-132-10310-215-13-1194-426-167-332386-8571726-3502-13-4253011610-11-917342-122462-71392012400-18-3-24-863-1];>>b=[5;12;3;2;3;46;13;38;19;-21];>>LUresolve(a,b)z=101l=1.00000000000002.00001.0000000000001.000001.000000000000-2.00003.00001.0000000000-1.00001.00004.0000-0.46671.0000000002.00001.0000-5.0000-0.2667-1.64411.000000000-1.00003.0000-0.0667-0.08470.29691.00000004.0000-1.00006.0000-0.2667-3.76272.73033.18631.0000001.0000-4.000026.00001.26675.0000-0.0182-18.535610.75921.000000-7.000030.00004.600021.7288-11.1148-59.930629.21790.98561.0000u=1.0e+03*0.00400.0020-0.0030-0.00100.00200.0010000000.00200.00100.00500.01000.00700.00200.00300.00200.0020000.001000.00200-0.0030-0.00200.0010-0.00100000.01500.01300.03100.04000.05400.06300.06500000-0.00390.01550.03410.07030.09270.1261000000.04170.05180.17280.33610.56170000000.0062-0.0297-0.1214-0.26720000000-0.0840-0.1669-0.349500000000-0.9987-1.8562000000000-0.7229y=1.0e+03*0.00500.0020-0.00200.01200.01960.05940.0058-0.07180.84281.8498x=8.377841.271410.527830.7380-28.04387.3550-14.84783.72733.9121-2.5589结果分析和讨论LU分解所得到结果x=(8.3778,41.2714,10.5278,30.7380,-28.0438,7.3550,-14.8478,3.7273,3.9121,-2.5589)T与正确解含有很大误差。这是因为系数矩阵本身数值性质所决定,因为计算过程中并未进行选主元过程,所以由系数矩阵产生L和U矩阵就含有了很大数值问题。由L和U所求出x和y就会产生很大误差。这是由矩阵本身数值问题所引发,与算法本身无关,消除误差关键在于计算过程中需要进行选主元。完成题目体会与收获对于其她解线性方程组方法来看,LU分解是较为高效,了解并熟练利用LU分解对于学习数值计算与解线性方程组都有很大帮助。在进行分解过程中应注意矩阵数值问题与怎样选择主元问题,这么才能得到方程组正确解,不然将产生非常大误差。在进行分解时应该格外注意,因为系数矩阵存在很多数值问题。（2）设对称正定阵系数阵线方程组数学原理平方根法解n阶线性方程组Ax=bcholeskly方法也叫做平方根法,这里对系数矩阵A是有要求,需要A是对称正定矩阵,依据数值分析相关理论,假如A对称正定,那么系数矩阵就能够被分解为形式,其中L是下三角矩阵,将其代入Ax=b中,可得:进行以下分解:那么就可先计算y,再计算x,因为L是下三角矩阵,是上三角矩阵,这么计算比直接使用A计算简便,同时你应该也发觉了工作量就转移到了矩阵分解上面,那么对于对称正定矩阵A进行Cholesky分解,我再描述一下过程吧:假如你对原理很清楚那么这一段能够直接跳过。设,即其中第1步,由矩阵乘法,故求得通常,设矩阵L前k-1列元素已经求出第k步,由矩阵乘法得于是改善平方根法在平方根基础上,为了避免开方运算,所以用计算;其中,;得按行计算元素及对元素公式对于..计算出第行元素后,存放在第行相置,然后再计算第行元素,存放在第行.对角元素存放在对应位置.对称正定矩阵按分解和按分解计算量差不多,但分解不需要开放计算。求解,计算公式分别以下公式。程序设计1、平方根法function[x]=pfpf(A,b)%楚列斯基分解求解正定矩阵线性代数方程A=LL’先求LY=b再用L’X=Y即能够求出解X[n,n]=size(A);L(1,1)=sqrt(A(1,1));fork=2:nL(k,1)=A(k,1)/L(1,1);endfork=2:n-1L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));fori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));%解下三角方程组Ly=b对应递推公式以下,求出y矩阵y=zeros(n,1);%先生成方程组因变量位置,给定y初始值fork=1:nj=1:k-1;y(k)=(b(k)-L(k,j)*y(j))/L(k,k);end%解上三角方程组L’X=Y递推公式以下,可求出X矩阵x=zeros(n,1);U=L';%求上对角矩阵fork=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j))/U(k,k);end>>A=[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,0-4,-1,14,1,-8,-3,5,60,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,-4,4,11,1,-40,2,5,-3,-10,1,14,20,0,6,3,-3,-4,2,19];>>b=[0;-6;20;23;9;-22;-15;45];>>x=pfpf(A,b)x=121.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.01852、改善平方根法function[x]=improvecholesky(A,b,n)%用改善平方根法求解Ax=bL=zeros(n,n);%L为n*n矩阵D=diag(n,0);%D为n*n主对角矩阵S=L*D;fori=1:n%L主对角元素均为1L(i,i)=1;endfori=1:nforj=1:n%验证A是否为对称正定矩阵if(eig(A)<=0)|(A(i,j)~=A(j,i))%A特征值小于0或A非对称时,输出wrongdisp('wrong');break;endendendD(1,1)=A(1,1);%将A分解使得A=LDLTfori=2:nforj=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);%x,y为n*1阶矩阵x=zeros(n,1);fori=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);%经过LDy=b解得y值endfori=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));%经过LTx=y解得x值end>>A=[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,0-4,-1,14,1,-8,-3,5,60,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,-4,4,11,1,-40,2,5,-3,-10,1,14,20,0,6,3,-3,-4,2,19];>>b=[0;-6;20;23;9;-22;-15;45];>>n=8;>>x=improvecholesky(A,b,n)x=121.1481-140.112729.7515-60.152810.9120-26.79635.4259-2.0185结果分析和讨论平方根法和改善平方根法求解线性方程组解为x=（121.1481,-140.1127,29.7515,-60.1528,10.9120,-26.7963,5.4259,-2.0185）T。与正确解相比较也存在很大误差,即使系数矩阵对角元素都大于零, 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 上能够无须选择主元,但因为矩阵数值问题较大,不选主元结果就是产生很大误差,所以在求解过程中还是应该选择主元以此消除误差,提升精度。完成题目体会与收获对称正定矩阵平方根法及改善平方根法是现在处理这类问题最有效方法之一,合理利用话,能够产生很好求解效果。改善平方根法较平方根法,因为不用进行开方运算,所以含有一定求解优势。经过求解此题,使我学会了怎样利用Matlab编程来处理平方根法和改善平方根法问题。（3）三对角形线性方程组数学原理设系数矩阵为三对角矩阵则方程组Ax=f称为三对角方程组。设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记　　　可先依次求出L,U中元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。实际上,求解三对角方程组2追赶法将矩阵三角分解计算与求解两个三角方程组计算放在一起,使算法更为紧凑。其计算公式为:(*)三、程序设计functionx=chase(a,b,c,f)%求解线性方程组Ax=f,其中A是三对角阵%a是矩阵A下对角线元素a(1)=0%b是矩阵A对角线元素%c是矩阵A上对角线元素c(n)=0%f是方程组右端向量n=length(f);x=zeros(1,n);y=zeros(1,n);d=zeros(1,n);u=zeros(1,n);%预处理d(1)=b(1);fori=1:n-1u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i);end%追过程y(1)=f(1)/d(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/d(i);end%赶过程x(n)=y(n);fori=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);end>>a=[0,-1,-1,-1,-1,-1,-1,-1,-1,-1];>>b=[4,4,4,4,4,4,4,4,4,4];>>c=[-1,-1,-1,-1,-1,-1,-1,-1,-1,0];>>f=[7,5,-13,2,6,-12,14,-4,5,-5];>>x=chase(a,b,c,f)x=2.00001.0000-3.00000.00001.0000-2.00003.0000-0.00001.0000-1.0000结果分析和讨论追赶法求解结果为x=（2,1,-3,0,1,-2,3,0,1,-1）T。求解结果与正确解一样,这表明追赶法对于求解三对角方程组含有非常高精度,误差非常小。算法次数也较少,不选主元也能够有效算出正确结果,是一个计算量少而数值稳定方法。完成题目体会与收获经过本题求解,加深了对追赶法求解三对角方程组算法原理了解。利用追赶法Matlab编程,并学会了又一个求解特殊方程组方法。在计算量方面,追赶法有着巨大优势,所以在可能情况下应优先使用追赶法。加深了对数值计算教材知识了解,收获非常大。题目四:微分方程数值解在传染病传染理论中,一个比较初等微分方程可用于估计在任何时刻人群中受传染数量,只要做了合适简化。尤其,假定在一个固定人群中,全部些人含有一样机会被感染,且一感染就保持这种状态。假设表示在时刻易受感染人数量,表示感染他人人数。有理由假设感染他人人数改变速率与和乘积成正比,因为速率既取决于感染他人人数也取决于那个时刻易感染人数。假如人口多足以假定和是连续变量,则问题表示为,其中是常数,,即为人口总数。这个方程就变为仅包含:问题:一个地域假定,又假定时间用天来衡量,求30天结束时感染他人人口数量近似值数学原理Eular方法:一阶线性微分方程初值问题（1）方程离散化:差分和差商（2）经过初始值,依据递推公式（2）逐步算出就为显性Eular方法。隐形Eular方法:（3）公式（3）即为隐式Eular公式。程序设计functionE2=Euler_2(fun,x0,y0,xN,N)%向后Euler公式,其中,%fun为一阶微分方程函数%x0,y0为初始条件%xN为取值范围一个端点%h为区间步长%N为区间个数%x为xN组成向量%y为yN组成向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;forn=1:N%用迭代法求y(n+1)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));fork=1:3z1=y(n)+h*feval(fun,x(n+1),z0);ifabs(z1-z0)<1e-3break;endz0=z1;endy(n+1)=z1;endT=[x',y']functionz=f(x,y)z=(2*10^(-6))*(100000-y)*y;>>Euler_2('f',0,1000,29,29)T=1.0e+04*00.10000.00010.12460.00020.15510.00030.19290.00040.23960.00050.29720.00060.36800.00070.45480.00080.56050.00090.68860.00100.84290.00111.02710.00121.24500.00131.49990.00141.79430.00152.12950.00162.50490.00172.91820.00183.36470.00193.83770.00204.32870.00214.82810.00225.32600.00235.81280.00246.28000.00256.72080.00267.13000.00277.50450.00287.84280.00298.1449结果分析和讨论用隐式欧拉法求得第30天时感染他人人口数量近似值为y（29）=81449人。隐式欧拉法较之欧拉法误差更小,求取结果含有一定可参考价值。对于结果存在误差,只能经过改变算法实现。完成题目体会与收获求解初值问题,欧拉法是最简单数值解法。而隐式欧拉法更深入降低了欧拉法误差。隐式欧拉法原理依旧十分简单。在对数值结果正确度并非要求很高时可采取隐式欧拉法,这么既能快速地得到结果,又能得到适宜结果。本题处理对我数值计算学习有着很大帮助,使我对初值问题了解愈加深刻,并熟悉了初值问题求解方法和求解过程,。题目五:数值积分在光学物理中研究矩角位光绕射常常会用到Fresnel积分:对于结构一个正确到值表格供研究者查询（Romberg积分算法）。数学原理龙贝格算法是由递推算法得来。由梯形公式得出辛普生公式得出柯特斯公式最终得到龙贝格公式。在变步长过程中探讨梯形法计算规律。设将求积区间[a,b]分为n个等分,则一共有n+1个等分点,n。这里用表示复化梯形法求得积分值,其下标n表示等分数。先考察下一个字段[],其中点,在该子段上二分前后两个积分值显然有下列关系将这一关系式相关k从0到n-1累加求和,即可导出下列递推公式需要强调指出是,上式中代表二分前步长,而梯形法算法简单,但精度低,收敛速度缓慢,怎样提升收敛速度以节省计算量,自然式大家极为关心。依据梯形法误差公式,积分值截断误差大致与成正比,所以步长减半后误差将减至四分之一,现有将上式移项整理,知由此可见,只要二分前后两个积分值和相当靠近,就能够确保计算确保结果计算结果误差很小,这种直接用计算结果来估量误差方法称作误差事后估量法。按上式,积分值误差大致等于,假如用这个误差值作为一个赔偿,能够期望,所得应该是愈加好结果。按上式,组合得到近似值直接验证,用梯形二分前后两个积分值和按式组合,结果得到辛普生法积分值。再考察辛普生法。其截断误差与成正比。所以,若将步长折半,则误差对应减至十六分之一。现有由此得不难验证,上式右端值其实就等于,就是说,用辛普生法二分前后两个积分值和,在按上式再做线性组合,结果得到柯特斯法积分值,现有反复一样手续,依据斯科特法误差公式可深入导出龙贝格公式应该注意龙贝格公式已经不属于牛顿—柯特斯公式范围。在步长二分过程中利用公式加工三次,就能将粗糙积分值逐步加工成精度较高龙贝格,或者说,将收敛缓慢梯形值序列加工成熟练快速龙贝格值序列,这种加速方法称龙贝格算法。程序设计functionlbg(fx,a,t,k,e)%fx为要求积分函数%a为要求积分下限%t为要求积分上限%k为要求积分最大次数%e为要求积分结束精度T=zeros(k,k);%T为龙贝格算法递推表T(1,1)=(t-a)*(1+fx(t))/2;fori=1:km=0;forj=1:2^(i-1)m=m+fx(a+(2*j-1)*(t-a)/(2^i));endT(i+1,1)=0.5*T(i,1)+(t-a)*m/2^i;forn=1:iT(i+1,n+1)=(4^n*T(i+1,n)-T(i,n))/(4^n-1);endifabs(T(i+1,i+1)-T(i,i))<=e&i>=4break;else;endendfori=1:kifT(i,1)==0j=i;break;else;endendifj==kerror('所求次数不够或不可积')else;endT=T(1:j-1,1:j-1)disp('所求积分值为:')%输出求得积分值disp(T(j-1,1))functionfx=f(x)fx=cos(pi*x^2/2);functionfx=f(x)fx=sin(pi*x^2/2);>>a=0;t=0.1;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.100000000.10000.10000000.10000.10000.1000000.10000.10000.10000.100000.10000.10000.10000.10000.1000所求积分值为:0.1000>>a=0;t=0.2;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.199800000.19990.19990000.19990.19990.1999000.19990.19990.19990.199900.19990.19990.19990.19990.1999所求积分值为:0.1999>>a=0;t=0.3;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.298500000.29920.29940000.29930.29940.2994000.29940.29940.29940.299400.29940.29940.29940.29940.2994所求积分值为:0.2994>>a=0;t=0.4;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.393700000.39650.39740000.39720.39750.3975000.39740.39750.39750.397500.39750.39750.39750.39750.3975所求积分值为:0.3975>>a=0;t=0.5;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.481000000.48930.49210000.49160.49230.4923000.49210.49230.49230.492300.49230.49230.49230.49230.4923所求积分值为:0.4923>>a=0;t=0.6;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.553300000.57370.58040000.57920.58110.5811000.58060.58110.58110.581100.58100.58110.58110.58110.5811所求积分值为:0.5810>>a=0;t=0.7;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.601300000.64420.65850000.65580.65960.6597000.65870.65960.65970.659700.65940.65970.65970.65970.6597所求积分值为:0.6594>>a=0;t=0.8;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.614300000.69460.72140000.71580.72280.7229000.72110.72280.72280.722800.72240.72280.72280.72280.7228所求积分值为:0.7224>>a=0;t=0.9;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.582300000.71860.76400000.75340.76500.7650000.76200.76480.76480.764800.76410.76480.76480.76480.7648所求积分值为:0.7641>>a=0;t=1;k=100;e=10^(-4);fx=@(x)cos(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.500000000.71190.78260000.76340.78050.7804000.77580.77990.77990.779900.77890.77990.77990.77990.7799所求积分值为:0.7789>>a=0;t=0.1;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.05080000000000.02560.0172000000000.01300.00890.008300000000.00680.00470.00440.00440000000.00360.00260.00250.00240.0024000000.00210.00160.00150.00150.00150.001500000.00130.00100.00100.00100.00100.00100.00100000.00090.00080.00080.00080.00080.00080.00080.0008000.00070.00070.00060.00060.00060.00060.00060.00060.000600.00060.00060.00060.00060.00060.00060.00060.00060.00060.0006所求积分值为:6.2125e-04>>a=0;t=0.2;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.106300000000000.05470.03750000000000.02930.02090.0197000000000.01670.01250.01200.011800000000.01040.00840.00810.00800.00800000000.00730.00630.00610.00610.00610.0061000000.00580.00520.00520.00510.00510.00510.005100000.00500.00470.00470.00470.00470.00470.00470.00470000.00460.00440.00440.00440.00440.00440.00440.00440.0044000.00440.00430.00430.00430.00430.00430.00430.00430.00430.004300.00430.00430.00420.00420.00420.00420.00420.00420.00420.00420.0042所求积分值为:0.0043>>a=0;t=0.3;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.171100000000000.09090.06410000000000.05210.03910.0375000000000.03300.02660.02580.025600000000.02350.02040.02000.01990.01980000000.01880.01720.01700.01700.01700.0170000000.01650.01570.01560.01560.01550.01550.015500000.01530.01490.01480.01480.01480.01480.01480.01480000.01470.01450.01450.01450.01450.01450.01450.01450.0145000.01440.01430.01430.01430.01430.01430.01430.01430.01430.014300.01430.01420.01420.01420.01420.01420.01420.01420.01420.01420.0142所求积分值为:0.0143>>a=0;t=0.4;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.2497000000000000.13740.100000000000000.08440.06670.06450000000000.05860.05000.04890.0487000000000.04590.04170.04110.04100.041000000000.03960.03750.03720.03720.03720.03720000000.03650.03540.03530.03530.03530.03530.0353000000.03490.03440.03430.03430.03430.03430.03430.034300000.03410.03390.03380.03380.03380.03380.03380.03380.03380000.03380.03360.03360.03360.03360.03360.03360.03360.03360.0336000.03360.03350.03350.03350.03350.03350.03350.03350.03350.03350.033500.03350.03340.03340.03340.03340.03340.03340.03340.03340.03340.03340.0334所求积分值为:0.0335>>a=0;t=0.5;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.3457000000000000.19730.147900000000000.12910.10640.10360000000000.09650.08560.08420.0839000000000.08050.07510.07450.07430.074300000000.07260.06990.06960.06950.06950.06950000000.06860.06730.06720.06710.06710.06710.0671000000.06670.06600.06590.06590.06590.06590.06590.065900000.06570.06540.06530.06530.06530.06530.06530.06530.06530000.06520.06510.06500.06500.06500.06500.06500.06500.06500.0650000.06500.06490.06490.06490.06490.06490.06490.06490.06490.06490.064900.06490.06480.06480.06480.06480.06480.06480.06480.06480.06480.06480.0648所求积分值为:0.0649>>a=0;t=0.6;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.4607000000000000.27260.209900000000000.18850.16050.15720000000000.14880.13550.13390.1335000000000.12950.12300.12220.12200.122000000000.12000.11680.11640.11630.11630.11630000000.11520.11370.11350.11340.11340.11340.1134000000.11290.11210.11200.11200.11200.11200.11200.112000000.11170.11130.11130.11130.11130.11130.11130.11130.11130000.11110.11090.11090.11090.11090.11090.11090.11090.11090.1109000.11080.11070.11070.11070.11070.11070.11070.11070.11070.11070.110700.11070.11060.11060.11060.11060.11060.11060.11060.11060.11060.11060.1106所求积分值为:0.1107>>a=0;t=0.7;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.59360000000000000.36370.2871000000000000.26370.23040.226600000000000.21690.0.19940.19890000000000.19430.18670.18570.18550.1855000000000.18310.17940.17890.17880.17880.178800000000.17760.17580.17550.17550.17550.17550.17550000000.17490.17400.17380.17380.17380.17380.17380.1738000000.17350.17300.17300.17300.17300.17300.17300.17300.173000000.17280.17260.17260.17260.17260.17260.17260.17260.17260.17260000.17250.17240.17230.17230.17230.17230.17230.17230.17230.17230.1723000.17230.17230.17220.17220.17220.17220.17220.17220.17220.17220.17220.172200.17220.17220.17220.17220.17220.17220.17220.17220.17220.17220.17220.17220.1722所求积分值为:0.1722>>a=0;t=0.8;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.73770000000000000.46830.3785000000000000.35390.31570.311600000000000.30050.28270.28050.28000000000000.27460.26600.26490.26470.2646000000000.26190.25770.25710.25700.25700.257000000000.25560.25350.25320.25320.25320.25310.25310000000.25250.25140.25130.25130.25120.25120.25120.2512000000.25090.25040.25030.25030.25030.25030.25030.25030.250300000.25010.24990.24980.24980.24980.24980.24980.24980.24980.24980000.24970.24960.24960.24960.24960.24960.24960.24960.24960.24960.2496000.24950.24950.24950.24950.24950.24950.24950.24950.24950.24950.24950.249500.24940.24940.24940.24940.24940.24940.24940.24940.24940.24940.24940.24940.2494所求积分值为:0.2494>>a=0;t=0.9;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=0.88010000000000000.58080.4810000000000000.45590.41430.409800000000000.39690.37720.37480.37420000000000.36810.35850.35730.35700.3569000000000.35390.34920.34850.34840.34840.348300000000.34680.34450.34420.34410.34410.34410.34410000000.34330.34210.34200.34190.34190.34190.34190.3419000000.34150.34090.34090.34090.34080.34080.34080.34080.340800000.34070.34040.34030.34030.34030.34030.34030.34030.34030.34030000.34020.34010.34000.34000.34000.34000.34000.34000.34000.34000.3400000.34000.33990.33990.33990.33990.33990.33990.33990.33990.33990.33990.339900.33990.33980.33980.33980.33980.33980.33980.33980.33980.33980.33980.33980.3398所求积分值为:0.3399>>a=0;t=1;k=100;e=10^(-4);fx=@(x)sin(pi*x^2/2);>>lbg(fx,a,t,k,e)T=1.00000000000000000.69130.5885000000000000.56340.52080.516300000000000.50080.47990.47720.47650000000000.46950.45910.45770.45740.4573000000000.45390.44870.44800.44780.44780.447800000000.44610.44350.44310.44300.44300.44300.44300000000.44220.44090.44070.44070.44060.44060.44060.4406000000.44020.43960.43950.43950.43950.43940.43940.43940.439400000.43920.43890.43890.43890.43890.43890.43890.43890.43890.43890000.43870.43860.43860.43860.43860.43860.43860.43860.43860.43860.4386000.43850.43840.43840.43840.43840.43840.43840.43840.43840.43840.43840.438400.43840.43830.43830.43830.43830.43830.43830.43830.43830.43830.43830.43830.4383所求积分值为:0.4384结果分析和讨论对于不一样t值求得不一样c(t)和s(t)值以下表所表示。t0.10.10006.2125e-040.20.19990.00430.30.29940.01430.40.39750.03350.50.49230.06490.60.58100.11070.70.65940.17220.80.72240.24940.90.76410.33991.00.77890.4384伴随t增大,c(t)和s(t)值也逐步变大,这符合积分改变原理,所得计算值能够被研究者所使用。完成题目体会与收获经过对此题编程计算,基础学会了使用龙贝格算法利用Matlab处理相关问题,同时加深了对于龙贝格计算原理了解。同时了解到在光学物理中研究矩角位光绕射常常会用到Fresnel积分。经过这次作业练习,使得我对数值计算相关知识概念有了不一样了解。附录牛顿插值程序function[y,A,C,L]=newdscg(X,Y,x,M)%y为对应x值,A为差商表,C为多项式系数,L为多项式%X为给定节点,Y为节点值,x为待求节点n=length(X);m=length(x);%n为X长度fort=1:mz=x(t);A=zeros(n,n);A(:,1)=Y';s=0.0;p=1.0;q1=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));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));fork=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)=polyval(C,z);%输出y值endL(k,:)=poly2sym(C);%输出多项式最小二乘法拟合程序>>x=[0,5,10,15,20,25,30,35,40,45,50,55];>>y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.0