数值分析实验题

数值实验题解题报告力0王天宝20100116321数值积分试用不同方法计算定积分I≡∫0π2xsin(2x)2cos2(x)+sin2(x)dx(1)使其误差≤10−7(此积分的准确值为πln√2+12)。(1)用复合Simpson求积公式，打印出子区间个数，以及定积分的近似值。(2)采用n=4的Gauss-Legendre求积公式的复合求积公式，打印出计算函数值的节点，以及定积分的近似值。由题意，积分函数f(x)=xsin(2x)1+cos2(x)。将区间[0,π/2]分为n等分，xk=kh,h=π2n,k=0,1,⋯,n.I)用复合Simpson公式求解在每个子区间[xk−1,xk]上采用Simpson求积公式，那么有I=∑k=1n∫xk−1xkf(x)dx=h6∑k=1n[f(xk−1)+4f(xk−1/2)+f(xk)]+̃En(f).其中xk−1/2=12(xk−1+xk).令Sn(f)=h6∑k=1n[f(xk−1)+4f(xk−1/2)+f(xk)]也可以写成Sn(f)=h6[f(0)+∑k=1n(4f(xk−1/2)+2f(xk))−f(π2)]根据上式用Matlab编程求解积分近似值。程序输出：子区间个数：50定积分近似值：0.591330706378764误差：1.063368537668197e-08II)用n=4的Gauss-Legendre求积公式的复合求积公式求解对于每个子区间[xk−1,xk]，令t=2h(x−xk−1)−1,则x=h2(t+1)+xk−1,dx=h2dt，那么I=∑k=1n∫xk−1xkf(x)dx=∑k=1n∫−11f(h2(t+1)+xk−1)h2dt再根据n=4的Gauss-Legendre积分公式，I=∑k=1n∑m=04Amf(h2(tm+1)+xk−1)h2+̃En(f)令Sn(f)=h2∑k=1n∑m=04Amf(h2(tm+1)+xk−1)(1.1)其中Gauss求积点t0,t1,⋯,t4为5次Legendre多项式P5的零点，P5(x)=18(63x5−70x3+15x)因此t0,t1,⋯,t4可根据方程63x5−70x3+15x=0(1.2)求解得到，相应的系数A0,A1,⋯,A4可由式Am=251P4(tm)P5'(tm),m=0,1,⋯,4(1.3)计算得到，其中P4(x)=18(35x4−30x2+3)P5'(x)=18(315x4−210x2+15)综上，用Matlab编程求解的流程为，先求解方程(1.2)，得到t0,t1,⋯,t4，再根据式(1.3)求出A0,A1,⋯,A4，最后用求和式(1.1)计算积分近似值Sn(f)，打印出积分节点A0,A1,⋯,A4。程序输出：子区间个数：50定积分近似值：0.591330695745079积分节点：0.0000000000000000.9061798459386640.538469310105683-0.906179845938664-0.538469310105683误差：5.551115123125783e-16比较二者的结果可以看出，子区间个数相同，n=4的Gauss-Legendre复合求积公式精度比Simpson复合求积公式的精度高得多。2病态线性方程组的求解考虑下面的线性代数方程组的解Hx=b,(2)其中H=(hij)n×n,hij=1i+j−1,1≤i,j≤n.b=(b1,b2,⋯,bn)T,bi=∑j=1nhij,1≤i≤n.方程组存在唯一解x∗=(1,1,⋯,1)T.(1)选取n=6，分别用Gauss消去法，J迭代方法，SOR迭代方法(ω=1,1.25,1.5等)求解方程组(2)，将结果与标准解比较。(2)逐步增大维数n(n=8,10,⋯)，重复求解方程组，将计算结果与准确解进行比较。结果说明什么？容易证明，J迭代法不收敛，因此只考虑Gauss消去法和SOR迭代法。Gauss消去法对于方程组(2)，需要对A=[H∣b]进行顺序消去。顺序消去第1步，若a11≠0，令li1=ai1a11,i=2,3,⋯,n.做运算aij=aij−li1a1j,i,j=2,3,⋯,n.得到新的A；若a11=0，则往下寻找am1≠0，交换A第1行和第m行，得到新的A满足a11≠0，再做上面的运算。顺序消去第k步，若akk≠0，令lik=aikakk,i=k+1,⋯,n.做运算aij=aij−likakj,i,j=k+1,⋯,n.得到新的A；若akk=0，则往下寻找am1≠0，交换A第k行和第m行，得到新的A满足akk≠0，再做上面的运算。回代过程，xn=an,n+1ann,xi=1aii(ai,n+1−∑j=i+1naijxj),i=n−1,⋯,2,1根据上述步骤可以编写Gauss消去法求解程序。SOR迭代法对于方程组(2)，记H=D−L−U，其中D=diag(h11,h22,⋯,hnn)是h的对角线部分，而L=−[000⋯0h2100⋯0h31h320⋯0⋮⋮⋱⋱⋮hn1hn2⋯hn,n−10],U=−[0h12h13⋯h1n00h23⋯h2n000⋱⋮⋮⋮⋮⋱hn−1,n000⋯0]SOR迭代法的向量形式为x(k+1)=(D−ωL)−1[(1−ω)D+ωU]x(k)+ω(D−ωL)−1b令M=1ω(D−ωL),N=1ω[(1−ω)D+ωU]则迭代方法为x(k+1)=M−1Nx(k)+M−1b根据上式编写SOR迭代法求解程序。I)程序输出：n=6时Gauss消去法求得方程的解为0.9999999999992281.0000000000219370.9999999998517921.0000000003853700.9999999995745841.000000000167680误差|x-x*|=6.16e-10n=6,松弛因子=1.00,误差要求|x-x*|=1.00e-02,SOR迭代法求得方程的解为0.9998858968532541.0016168852599310.9957475252246670.9993638527280981.0099998383629900.993293213357308迭代次数：14508n=6,松弛因子=1.25,误差要求|x-x*|=1.00e-02,SOR迭代法求得方程的解为1.0000194436607160.9994610256969581.0035742412884580.9908381664641551.0099999903701220.996094540492699迭代次数：82840n=6,松弛因子=1.50,误差要求|x-x*|=1.00e-02,SOR迭代法求得方程的解为1.0000203026437480.9994465015869691.0036285353156680.9907777978411571.0099999965222600.996114916576111迭代次数：174587n=8,Gauss消去法误差|x-x*|=7.09e-07n=10,Gauss消去法误差|x-x*|=7.08e-04n=12,Gauss消去法误差|x-x*|=3.31e-01n=14,Gauss消去法误差|x-x*|=1.71e+01n=16,Gauss消去法误差|x-x*|=1.57e+01n=18,Gauss消去法误差|x-x*|=6.90e+01n=20,Gauss消去法误差|x-x*|=1.05e+02松弛因子=1.00,误差要求|x-x*|=1.00e-02,对于不同的n,SOR迭代法迭代次数统计n=8,迭代次数：42694n=10,迭代次数：25508n=12,迭代次数：88174n=14,迭代次数：59866n=16,迭代次数：46548n=18,迭代次数：44716n=20,迭代次数：144587n=100,迭代次数：118156根据上面的输出可以得出结论：Gauss消去法对于病态线性方程组的求解结果会随着维数的增大越来越发散，n比较大时所求的结果与准确解相差甚远；SOR迭代法求解病态线性方程组时，只要迭代次数足够，可以得到足够精确的解，一般来说n越大，迭代次数越多，在n=6时，松弛因子取为1.00时迭代收敛最快。3非线性方程组的求解已知非线性方程组{3x1−cos(x1x2)−12=0,x13−81(x2+0.1)2+sin(x3)+1.06=0,e−x1x2+20x3+10π−33=0.(3)(1)建立一个在区域D={(x1,x2,x3)∣∣xi∣≤1,i=1,2,3}上满足压缩映射定理的不动点迭代法，计算方程组的解。(2)用Newton法求解方程组(3)，试验不同的初值。I)不动点迭代法将方程组化为等价的x=Φ(x),其中x=[x1x2x3],Φ(x)=[ϕ1(x)ϕ2(x)ϕ3(x)]=[13(cos(x1x2)+12)±19√1.06+sin(x3)+x13−0.1120(−e−x1x2−10π−33)]很容易验证Φ在D上的压缩性。分别令ϕ2(x)=19√1.06+sin(x3)+x13−0.1和ϕ2(x)=−19√1.06+sin(x3)+x13−0.1，取x0=(0,0,0)T,根据迭代计算公式xk+1=Φ(xk)进行迭代至误差ε=∥xk+1−x∥≤2.22×10−16.程序输出：(x1和x2为两个解，n1和n2为对应的迭代次数)x1=0.5000-0.0081-0.5238n1=12x2=0.4985-0.1916-0.5286n2=13II)Newton法将方程组化为等价的F(x)=0,其中F(x)=[3x1−cos(x1x2)−12x13−81(x2+0.1)2+sin(x3)+1.06e−x1x2+20x3+10π−33]从而F'(x)=[3+x2sin(x1x2)x1sin(x1x2)03x12−162x2−16.2cos(x3)−x2e−x1x2−x1e−x1x220]Newton迭代计算公式为xk+1=xk−(F'(xk))−1F(xk)迭代至ε=∥xk+1−xk∥≤2.22×10−16.取迭代初值x0=(0,0,0)T,程序输出：(x为解，n为迭代次数)x=0.5000-0.0081-0.5238n=6取迭代初值x0=(−1,−1,−1)T,程序输出：x=0.4985-0.1916-0.5286n=9取迭代初值x0=(1,1,1)T,x=0.5000-0.0081-0.5238n=9这样，也得到了两个解，与不动点迭代法求解结果相同。4常微分方程初值问题用四阶Runge-Kutta方法和Adams四阶预估-校正方法求解常微分方程{y'(x)=(x+2x3)y3(x)−xy(x),0x≤3y(0)=13(4)取步长h=0.02.将得到的数值结果在x=1,1.5,3处和问题的准确解y(x)=(3+2x2+6ex2)−12比较。需要求解的常微分方程为y'=f(x,y)其中f(x,y)=(x+2x3)y3−xy经典Runge-Kutta四阶方法{k1=f(xn,yn),k2=f(xn+h2,yn+h2k1),k3=f(xn+h2,yn+h2k2),k4=f(xn+h,yn+hk3),yn+1=yn+16h(k1+2k2+2k3+k4).(4.1)由初值条件，令y0=13，根据(4.1)式进行编程。Adams四阶预估-校正方法{yn+4[0]=yn+3+h24[55f(xn+3,yn+3)−59f(xn+2,yn+2)+37f(xn+1,yn+1)−9f(xn,yn)],fn+4[0]=f(xn+4,yn+4[0]),yn+4=yn+3+h24[9fn+4[0]+19f(xn+3,yn+3)−5f(xn+2,yn+2)+f(xn+1,yn+1)].(4.2)令y0=13，y1,y2,y3由runge-Kutta四阶方法求出，在根据(4.2)式，可以求出数值结果。程序输出：在x=1处，ans=0.21662642080.21662641710.2166264208在x=1.5处，ans=0.12458564