数值分析上机作业
数值分析上机作业(总)
数值分析上机实验一、解线性方程组直接法（教材49页14题）追赶法程序如下：function x=followup(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end主程序如下：function zhunganfaA=[2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0 0 -2 5 -2 0 0 0;0 0 0 -2 5 -2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5];b=[220/27;0;0;0;0;0;0;0];x=followup(A,b)计算结果：x =8.14784.07372.03651.01750.50730.25060.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor i=1:1:nA(i,j)=(1+0.1*i)^(j-1);endendc=cond(A)d=rcond(A)当n=5时c =5.3615e+005d =9.4327e-007当n=10时c =8.6823e+011d =5.0894e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。
数值分析大作业
数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
数值分析报告上机题课后作业全部-东南大学
实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。
按从小到大的顺序，有效位数分别为：5，6，6。
可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。
因此，采取从小到大的顺序累加得到的结果更加精确。
2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。
（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。
试确定尽可能大的δ。
○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。
算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。
运行Find.m，得到在不同的搜索精度下的最大sigma值。
数值分析上机作业
第一章第二题(1) 截断误差为104-时：k=1;n=0;m=0;x=0;e=1e-4;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141792613595791m =5001（2）截断误差为108-时：k=1;n=0;m=0;x=0;e=1e-8;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141592673590250m =50000001由以上计算可知，截断误差小于104-时，应取5001项求和，π=3.141792613595791；截断误差小于108-时，应取50000001项求和，π=3.141592673590250。
第二章第二题a=[0 -2 -2 -2 -2 -2 -2 -2];b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0];v=220;r=27;d=[v/r 0 0 0 0 0 0 0];n=8;for i=2:na(i)=a(i)/b(i-1);b(i)=b(i)-c(n-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(n)=d(n)/b(n);for i=n-1:-1:1d(i)=(d(i)-c(i)*d(i+1));end;I=d'I =1.0e+002 *1.490717294184090.704617906351300.311568212434910.128623612390290.049496991380330.017168822994210.004772412363470.00047741598598第三章第一题(1)Jacobi迭代法：b=[12;-27;14;-17;12]x = [0;0;0;0;0;]k = 0;r = 1;e=0.000001A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A));B = inv(D)*(D-A);f = inv(D)*b;p = max(abs(eig(B)));if p >= 1'迭代法不收敛'returnendwhile r >ex0 = x;x = B*x0 + f;k = k + 1;r = norm (x-x0,inf);endxk计算结果：x =1.0000-2.00003.0000-2.00001.0000k =65(2) 高斯赛德尔迭代：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;]x=[0;0;0;0;0];b=[12;-27;14;-17;12]c=0.000001L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-L)*U*x+inv(D-L)*b;k=1;while norm(X-x,inf)>= cx=X;X=inv(D-L)*U*x+inv(D-L)*b;k=k+1;endXk计算结果：X =1.0000-2.00003.0000-2.00001.0000k =37(3) SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0];b=[12;-27;14;-17;12]e=0.000001w=1.44;L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*bn=1;while norm(X-x,inf)>=ex=X;X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b;n=n+1;endXn计算结果：X =1.0000-2.00003.0000-2.00001.0000n =22由以上可知，共轭梯度法收敛速度明显快于Jacobi法和G-S法。
数值分析上机作业
第一章（1）按从大到小的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=2;j<=n;j++)s=s+1/(j*j-1);printf("s=%f\n",s);}(2)按从小到大的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=n;j>=2;j--)s=s+1/(j*j-1);printf("s=%f\n",s);}(3)精确值，从大到小，从小到大的计算程序为：#include<stdio.h>void main(){float n,j,s1,s2=0,s3=0;printf("请输入n：");scanf("%f",&n);s1=(3.0/2-1/n-1/(n+1))/2;for(j=2;j<=n;j++)s2=s2+1/(j*j-1);for(j=n;j>=2;j--)s3=s3+1/(j*j-1);printf("精确算法的结果为：%f\n",s1);printf("从大到小顺序的结果为：%f\n",s2); printf("从小到大顺序的结果为：%f\n",s3); }结果如下：精确值从大到小的值从小到大的值有效位数从大到小从小到大210S0.740049 0.740049 0.74005 6 5410S0.7499 0.749852 0.7499 4 4610S0.749999 0.749852 0.749999 3 6(4)通过本上机题，我们可以看出：当n较小时，两种算法的结果都很接近精确值，而当n较大时，两种算法的结果差别较大，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。
数值分析上机第四次作业
数值分析上机第四次作业实验项目：共轭梯度法求解对称正定的线性方程组实验内容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵，A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,nb[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
（3）*考虑模型问题，方程为222222(),(,)(0,1)(0,1)(0,)1,(1,),01(,0)1,(,1),01xy y x u u x y e x y D x yu y u y e y u x u x e x ∂∂+=+∈=⨯∂∂==≤≤==≤≤用正方形网格离散化，若取1/,10h N N ==，得到100n =的线性方程组，并用共轭梯度法（CG 法）求解，并对解作图。
实验要求：迭代初值可以取01(,1,...,)ij u i j N ==，计算到32||||10k r -≤停止．本题有精确解(,)xy u x y e =，这里k u 表示以k ij u 为分量的向量，u 表示在相应点(,)i j 上取值作为分量的向量．实验一：（1）编制函数子程序CGmethod 。
function [x,k]=CGmethod(A,b)n=length(A);x=zeros(n,1);r=b-A*x;rho=r'*r;k=0;while rho>10^(-12) & k<20k=k+1;if k==1p=r;elsebeta=rho/rho1;p=r+beta*p;endw=A*p;alpha=rho/(p'*w);x=x+alpha*p;r=r-alpha*w;rho1=rho;rho=r'*r;end编制主程序shiyan1_1：clear,clcA=[2,-1,0,0;-1,3,-1,0;0,-1,4,1;0,0,-1,5]; b=[3,-2,1,5]';[x,k]=CGmethod(A,b)运行结果为：x =1.3882-0.2855-0.02220.9367k =20（2）编制函数子程序CGmethod_1function [x,k]=CGmethod_1(A,b)n=length(A);x(1:n,1)=0;r=b-A*x;r1=r;k=0;while norm(r1,1)>=10^(-7)&k<10^4k=k+1;if k==1p=r;elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p;endr=r1;w=A*p;alpha=(r'*r)/(p'*w);x=x+alpha*p;r1=r-alpha*w;end编制主程序shiyan1_2:clear,clcn=1000;A=hilb(n);b=sum(A')';[x,k]=CGmethod_1(A,b)运行结果为：x 的值,均接近1，迭代次数k=32实验二实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
高等数值分析上机作业
画图分析
4
第 10 章 常微分方程的初边值问题解法
上机作业 2：
1. 分别用 Euler 方法（h=0.025）、改进的 Euler 方法（h=0.05）、经
典的 R-K 方法（h=0.1），求初值问题的数值计算 0.1,0.2,0.3,0.4,0.5
处解的近似值。
y y x2 1 y(0) 0.5
主程序： >> fun=inline('y-x.^2+1'); >> [x,y]=Euler1(fun,0.5,0.1)
运行结果：
x=
0.1000 0.2000 0.3000 0.4000 0.5000
y=
0.5000 0.6540 0.8200 0.9970 1.1841
用经典的R-K方法（）求解一阶微分方程组的初值问题，并分别画出
3
h=0.01;x=-1:h:1; B(1,1)=quad('asin(x)./(sqrt(1-x.^2))',-1,1); B(2,1)=quad('asin(x).*x./sqrt(1-x.^2)',-1,1); B(3,1)=quad('asin(x).*(2*x.^2-1)./(sqrt(1-x.^2))',-1,1); B(4,1)=quad('asin(x).*(4*x.^3-3*x)./(sqrt(1-x.^2))',-1,1); B A=diag([pi,pi/2,pi/2,pi/2]); c=inv(A)*B y=asin(x); plot(x,y,'r') hold on z=0.9204*x+0.4296*x.^3; plot(x,z,'k') hold on s=1.2733*x+0.1415*(4*x.^3-3*x); plot(x,s,'c') legend('y=arcsin(x)','y=0.9204x+0.4296x^3','s=1.2733x+0.141 5(4x^3-3*)') 最后可以得到切比雪夫拟合多项式为 y 1.2733x 0.1415(4x3 3x) ，
数值分析上机实验6
数值分析上机实验6根据表中数据，预测公元2000年时的世界人口。
问题分析与数学模型设人口总数为N(t)，根据人口理论的马尔萨斯模型，采用指数函数N(t) = e a + b t=+，令对数据进行拟合。
为了计算方便，将上式两边同取对数，得ln N a bty = ln N或N = e y变换后的拟合函数为y(t) = a + b t根据表中数据及等式 k k( 1，2，……，9)可列出关于两个未知数、b的9个方程的超定方程组（方程数多于未知数个数的方程组）a + t j b = y j(j= 1，2， (9)可用最小二乘法求解。
算法与数学模型求解算法如下：第一步：输入人口数据，并计算所有人口数据的对数值；第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；第三步：求解超定方程组并输出结果：a，b；第四步：利用数据结果构造指数函数计算2000年人口近似值N(2000)，结束。
MATLAB程序t=1960:1968;t0=2000;N=[29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.83];y=log(N);A=[ ones(9,1), t' ];d=A\ y' ;a=d(1),b=d(2)N0=exp(a+b*t0)x=1960:2001;yy=exp(a+b*x);plot(x,yy,t,N,'o',2000,N0,'o')计算结果为a =－3，b =6N (2000)所以取五位有效数，可得人口数据的指数拟合函数t e t N 0186.00383.33)(+-=经计算得2000年人口预测值为： (亿)。
例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier 级数的部分和（有限项）做拟合函数。
数值分析上机作业
数值分析上机作业2实验一：（1） ①用不动点迭代法求()013=--=x x x f 的根，至少设计两种迭代格式，一个收敛一个发散，1210-=ε。
（2） ②对迭代格式使用Aitken 加速，观察敛散性变化。
1取递推公式31)1(+=x x ，可以得到收敛时的迭代结果为：x=(2)^(1/3); t=1; while(1)if(abs(x-(x+1)^(1/3))<10^-12) break; endx=(x+1)^(1/3);t=t+1;end t xtemp=x^3-x-1 %带回来验证下 t = 16 x =1.324717957243755 temp =-4.225952920933196e-12 加速后代码如下 x=1;x=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)) t=1; while(1)if(abs(x-(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)))<10^-10) break; endx=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)); t=t+1; if (t>100000)break; %防止运算次数过多 end endfprintf('需要%d 次',t);输出需要670次>>此处取10^10-=ε，若到10^-12次方则可能需要运行更多次取13-=x x 则迭代发散。
使用aitken 加速计算结果如下 x=2; t=1; while(1)if( abs(x-((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1))<10^-10) break; end t=t+1;x=((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1); if (t>100000)break; %防止运算次数过多 end end t xt = 108x = 1.324717956244172由此可见经过aitken 加速以后，原来发散的迭代格式收敛了。
数值分析上机作业
数值分析上机作业（1、2、3、4、6章）第一章17. 舍入误差与有效数设2211NN j S j ==-∑，其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。
（1）编制按从大到小的顺序22211121311N S N =+++---，计算N S 的通用程序； （2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算N S 的通用程序；（3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数（编制程序时使用单精度）；（4）通过本上机题你明白了什么？运行结果：按从大到小的顺序计算得：N N S误差e有效位数2⨯8 100.7400495 94.95049501392230710-4⨯ 4 100.7498521 54.79049995000258010-6⨯ 3 100.7498521 41.46900000499994310-按从小到大的顺序计算得：N N S误差e有效位数2⨯84.95049501392230710-100.7400495 944.99950003618465610-⨯8 100.7499000 965.00044450291170510-⨯11 100.7499990 13（4）通过本题可以看出，不同算法造成的误差是不同的，好的算法可以让计算结果精度更高。
对于本题，当采用从大到小的顺序累加计算时，计算误差随着N的增大而增大；当采用从小到大的顺序累加计算时，计算误差随着N的增大而减小。
因此在N比较大时宜采用从小到达的顺序累加计算。
第二章20.Newton 迭代法（1）给定初值0x 及容许误差ε，编制Newton 法解方程()0f x =根的通用程序；（2）给定方程3()03x f x x =-=，易知其有三个根*1x =，*20x =，*3x =①由Newton 方法的局部收敛性可知存在0δ>，当0(,)x δδ∈-时Newton 迭代序列收敛于根*2x ，试确定尽可能大的δ；②试取若干初始值，观察当0(,1)x ∈-∞-，(1,)δ--，(,)δδ-，(,1)δ，(1,)+∞时Newton 序列是否收敛以及收敛于哪一个根；（3）通过本上机题，你明白了什么？本实验取610ε-=，找到的最大的0.774597δ=②当0(,1)x ∈-∞-时，计算结果如下：x0 xend -1000 -1.732051 -500 -1.732051 -100 -1.732051 -10 -1.732051 -1.5-1.732051Newton 序列收敛于-1.732051当0(1,)xδ∈--时，计算结果如下：x0 xend-0.95 1.732051-0.90 1.732051-0.85 1.732051-0.80 -1.732051-0.78 -1.732051 Newton序列收敛于1.732051或-1.732051当0(,)xδδ∈-时，计算结果如下：x0 xend-0.77 0.000000-0.5 0.000000-0.1 0.0000000.3 0.0000000.77 0.000000 Newton序列收敛于0当0(,1)xδ∈时，计算结果如下：x0 xend0.78 1.7320510.80 1.7320510.85 -1.7320510.90 -1.7320510.95 -1.732051 Newton序列收敛于1.732051或-1.732051当0(1,)x∈+∞时，计算结果如下：x0 xend1.5 1.73205110 1.732051100 1.732051500 1.7320511000 1.732051Newton序列收敛于1.732051（3）通过本题发现，Newton迭代法解方程初始值的选取非常重要，不同的初始值会收敛于方程不同的根，且有些区间是全局收敛，有些区间是局部收敛。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
昆明理工大学工科研究生《数值分析》上机实验
学院：材料科学与工程学院
专业：材料物理与化学
学号：2011230024
姓名: 郑录
任课教师：胡杰
P277-E1
1.已知矩阵A=
10787
7565
86109
75910
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦,B=
23456
44567
03678
00289
00010
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
,错误！未找到引用源。
=
11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
（1）用MA TLAB函数“eig”求矩阵全部特征值。
（2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。
解：MA TLAB程序如下：
求矩阵A的特征值：
clear;
A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
E=eig(A)
输出结果：
求矩阵B的特征值：
clear;
B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0];
E=eig(B)
输出结果：
求矩阵错误！未找到引用源。
的特征值：
clear; 错误！未找到引用源。
=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。
)
输出结果：
（2）A=
10
7877565861097
5
9
10
第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到：
第二部：[Q1,R1]=qr(A1);A2=R1*Q1
第三部：
[Q2,R2]=qr(A2);A3=R2*Q2
现在收缩，继续对A3的子矩阵错误！未找到引用源。
=29.8329
3.44110.00003.4411
4.30530.16110
0.1611
0.8516
-- 进行累世变换，得到（假设收缩后的矩阵为C6）
C6=29.8329
3.44110.00003.4411
4.30530.16110
0.1611
0.8516
-- 这是进行了6步qr 算法所得的结果。
故求的A 的近似特征值为错误！未找到引用源。
30.2886，错误！未找到引用源。
，错误！未找到引用源。
，错误！未找到引用源。
0.0102。
而A 的特
征值是错误！未找到引用源。
0.0102错误！未找到引用源。
30.2886
同理，用类似的方法可求矩阵B 和错误！未找到引用源。
的特征值，但过程过于繁琐，不再一一求解。
2用梯形公式、辛普森和Newton-Cotes 求积公式计算定积分⎰
π=
2
/0
I e x sin d x ，取精度
4
10
-，作出它们的积分图，并与精确值进行比较;
解 （1）用梯形求积公式计算定积分. 输入程序
>> h=pi/500; x=0:h:pi/2; y=exp(sin(x));
zt=trapz(x,y), ztc=cumtrapz(x,y), plot(x, ztc,'ro')
运行后屏幕显示用函数trapz 和cumtrapz 分别计算结果zt 、ztc 分别如下
zt =
3.10437572798742 ztc =
Columns 1 through 3
0 0.00630298652792 0.01264569951380
………………………………………………………………………..
Columns 250 through 251
3.08729642810745 3.10437572798742
（2）用辛普森求积公式计算定积分. 输入程序
>> syms x
L= inline(' exp(sin(x))');
[QS,FCNTS] =quad(L,0, pi/2,1.e-4,2)
运行后屏幕显示用辛普森求积公式计算定积分的值QS和递归次数FCNTS分别如下QS = FCNTS =
3.10438133817254 13
（3）用Newton-Cotes求积公式计算定积分. 在MATLAB6.5中输入程序
>> syms x
L= inline(' exp(sin(x))');
[Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3)
运行后屏幕显示用Newton-Cotes求积公式计算定积分的值Q8和递归次数FCNTS分别如下Q8 = FCNT8 =
3.10437901785555 33
（4）输入求定积分的精确值的程序
>> syms x
y=exp(sin(x)); F=int(y,0,pi/2), Fj=double(F)
wzt=abs( Fj- zt), wQS = abs(Fj- QS), wQ8 = abs(Fj- Q8)
运行后屏幕显示计算的定积分值F和近似值F j，梯形公式、辛普森和Newton-Cotes求积公式计算定积分的值与F j的绝对误差wztc, wQS和wQ8如下
Warning: Explicit integral could not be found.
> In C:\MATLAB6p5p1\toolbox\symbolic\@sym\int.m at line 58
F =
int(exp(sin(x)),x = 0 .. 1/2*pi)
Fj = wzt =
3.10437901785556 3.289868133471430e-006
wQS = wQ8 =
2.320316987436399e-006 7.993605777301127e-015
（5）输入画图程序
>> syms x
F=int(exp(sin(x)),0,pi/2),Fj=double(F),plot(Fj,'ro'),
hold on
L= inline(' exp(sin(x))');
Q=quad(L,0, pi/2,1.e-4,2),plot(Q,'g*')
hold off,hold on, h=pi/40;
x=0:h:pi/2; y=exp(sin(x));
ztc=cumtrapz(x,y),
plot(x, ztc,'mp'), hold off
hold on, [Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3), hold off
grid, xlabel('自变量 X'), ylabel('因变量 Y')
title('Newton-Cotes公式和梯形公式计算数值积分')
legend('精确值', ' 辛普森公式计算定积分','梯形公式计算定积分','Newton-Cotes公式计算定积分')
运行后屏幕显示图形（略）.。