矩阵与数值分析Word下载.docx

1、 T（i,j）=2; end if（i=（j+1） T（i,j）=-1; if（i=（j-1）I=eye（n-1）;k=1;for mm=1:（n-1） A（k:（k+n-2）,k:（k+n-2）=T+2*I; k=k+n-1;for xx=1:（n-2）（k+n-2）,（k+n-1）:（k+2*n-3）=-I; A（k+n-1）:（k+2*n-3）,k:（k+n-2）=-I;b=randn（n-1）2,1）;x0=zeros（n-1）2,1）;Jacobi迭代法Matlab程序如下：function n=jacobi（A,b,x0）eps= 1.0e-10;M = 10000;D=diag（

2、diag（A）; %求A的对角矩阵L=-tril（A,-1）; %求A的下三角阵U=-triu（A,1）; %求A的上三角阵B=D（L+U）;f=Db;x=B*x0+f;n=1; %迭代次数while norm（x-x0）=eps x0=x; x=B*x0+f; n=n+1; if（n=M） disp（Warning: 迭代次数太多，可能不收敛！）; return;Gauss-Seidel迭代法Matlab程序如下：function n=gauseidel（A,b,x0）G=（D-L）U;f=（D-L）b;x=G*x0+f; x=G*x0+f;A,b,x0=jz（10）;J10=jacobi（

3、A,b,x0）G10=gauseidel（A,b,x0）A,b,x0=jz（20）;J20=jacobi（A,b,x0）G20=gauseidel（A,b,x0）A,b,x0=jz（30）;J30=jacobi（A,b,x0）G30=gauseidel（A,b,x0） J10 =436；G10 =214J20 =1810；G20 =913J30 =3990；G30 =2001N=10且M=1000时，Jacobi迭代法和Gaussseidel迭代法的迭代次数分别为436和214；N=20且M=1000时，Jacobi迭代法和Gaussseidel迭代法的迭代次数分别为1810和913；N=30

4、且M=1000时 ，Jacobi和Gauss-seidel算法的迭代次数分别为3990和2001次。从以上结果可知在该问题下Jacobi算法的收敛性没有Gauss-seidel算法的收敛性好，原因在于Jacobi算法迭代过程中，对已算出的信息未加充分利用，一般来说，后面计算的值要比前面的计算值要精确些，而Gauss-seidel算法则充分利用了已求出来的信息，所以此算法的收敛性更为好一些。 2. 用Gauss列主元消去法、QR方法求解如下方程组：Gauss列主元消去法Matlab程序如下：function x,XA=GaussXQLineMain（A,b） N = size（A）;n = N（

5、1）;index = 0; me = max（abs（A（1:n,i）; %选取列主元 for k=i:n if（abs（A（k,i）=me） index = k; %保存列主元所在的行 break; temp = A（i,1:n）; A（i,1:n） = A（index,1: A（index,1:n） = temp; bb = b（index）; b（index）=b（i）; b（i） = bb; %交换主行 for j=（i+1）: if（A（i,i）=0）对角元素为0！ l = A（j,i）; m = A（i,i）; A（j,1:n）=A（j,1:n）-l*A（i,1:n）/m; b（j

6、）=b（j）-l*b（i）/m; %消元for i=n:1 if（i0）两端点函数值乘积大于0!else tol=1; fa=subs（sym（f）,findsym（sym（f）,a）; fb=subs（sym（f）,findsym（sym（f）,b）; root=a-（b-a）*fa/（fb-fa）;while（toleps） r1=root; fx=subs（sym（f）,findsym（sym（f）,r1）; s=fx*fa; if（s=0） root=r1; if（s root=b-（r1-b）*fb/（fx-fb）; root=a-（r1-a）*fa/（fx-fa）; tol=abs

7、（root-r1）;利用matlab的绘图功能可以确定求根区间为（1,3）root=Secant（exp（x）+2*x2+2*sin（x）-log（x）-16,1,3）root =1.962922683164462利用Newton迭代法求多项式的所有实零点，注意重根的问题。Newton迭代法的Matlab程序如下：x=-3;d=-1;while abs（d-x）0.5*10（-8） %在区间（-3，-1） d=x; x=x-（x4-3*x3-3*x2+11*x-6）/（4*x3-9*x2-6*x+11）;%牛顿法x0=dx=0;d=2;0.5*10（-8）%在区间（0，2）的重根 x=x-2*

8、（x4-3*x3-3*x2+11*x-6）/（4*x3-9*x2-6*x+11）;x1=dx=2.5;d=4;0.5*10（-8）%在区间（2.5，4）x2=dx0 =-2.00000000151233x1 =1.00000000065769x2 =3.00000000000002由以上结果可知，newton迭代法所得结果与初始值的选取密切相关。不同的初始值会产生不同的结果，甚至会影响牛顿迭代法收敛性。四、数值积分用三点Gauss型求积公式计算积分n点Gauss型求积公式matlab程序如下：function I = question4（f,a,b,n）%f为求积方程;%a,b分别为积分上下限

9、;%n为Gauss点个数;XK,AK=grule（n）;XK1=（b-a）/2）*XK+（a+b）/2）;I=（b-a）/2）*sum（AK.*subs（f,findsym（f）,XK1）;function bp,wf=grule（n）% bp,wf=grule（n）bp=zeros（n,1）; wf=bp; iter=2; m=fix（n+1）/2）; e1=n*（n+1）;mm=4*m-1; t=（pi/（4*n+2）*（3:4:mm）; nn=（1-（1-1/n）/（8*n*n）;xo=nn*cos（t）;for j=1:iter pkm1=1; pk=xo; for k=2: t1=x

10、o.*pk; pkp1=t1-pkm1-（t1-pkm1）/k+t1; pkm1=pk; pk=pkp1; den=1.-xo.*xo; d1=n*（pkm1-xo.*pk）; dpn=d1./den; d2pn=（2.*xo.*dpn-e1.*pk）./den; d3pn=（4*xo.*d2pn+（2-e1）.*dpn）./den; d4pn=（6*xo.*d3pn+（6-e1）.*d2pn）./den; u=pk./dpn; v=d2pn./dpn; h=-u.*（1+（.5*u）.*（v+u.*（v.*v-u.*d3pn./（3*dpn）; p=pk+h.*（dpn+（.5*h）.*（d

11、2pn+（h/3）.*（d3pn+.25*h.*d4pn）; dp=dpn+h.*（d2pn+（.5*h）.*（d3pn+h.*d4pn/3）; h=h-p./dp; xo=xo+h;bp=-xo-h;fx=d1-h.*e1.*（pk+（h/2）.*（dpn+（h/3）.*（. d2pn+（h/4）.*（d3pn+（.2*h）.*d4pn）;wf=2*（1-bp.2）./（fx.*fx）;if （m+m） n, bp（m）=0;if （m+m） = n）, m=m-1;jj=1:m; n1j=（n+1-jj）; bp（n1j）=-bp（jj）; wf（n1j）=wf（jj）;% endsyms

12、 x;f=exp（-x）;a=0;b=1;n=3;I = question4（f,a,b,n）I =0.632120255664068五、插值与逼近 1给定上的函数，请做如下的插值逼近：（1）构造等距节点分别取的Lagrange插值多项式；（2）取Chebyshev多项式的零点：作插值节点构造的插值多项式和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。（1）Lagrange插值matlab程序如下：function f = L1（n）syms t s;q=1/（1+s2）; x（i）=-1+2*（i-1）/（n-1）; y（i）=subs（q,s,x（i）;f = 0

13、.0;for（i = 1: l = y（i）; for（j = 1:i-1） l = l*（t-x（j）/（x（i）-x（j）; end; for（j = i+1: %计算拉格朗日基函数 f = f + l; %计算拉格朗日插值函数 f=collect（f）; %化简X=-1:0.1:1;Y=subs（q,s,X）;F=subs（f,t,X）;plot（X,Y,g,X,F,rf=collect（f）; %展开digits（5）;f=vpa（f）;1、在matlab命令窗口中输入：f = L1（5）运行结果2、在matlab命令窗口中输入：f = L1（8）3、在matlab命令窗口中输入：f

14、=L1（10）由图像可以看出n=5时拟合曲线与原始曲线的误差比较大，而n=8,10时拟合曲线与原始曲线几乎重合在一起。可见节点个越多，拉格朗日插值法所得的多项式曲线与原始曲线越接近，误差就越小。的零点为插值点,此时Lagrange插值matlab程序主体基本不变，变的只是插值点，其程序如下：function f = L2（n） x（i）=cos（2*i-1）*pi/2/n）; %计算拉格朗日插值函数 %绘图krpf = L2（10）有图像比较可以看出，N=10时，chebyshev零为拉格朗日插值时，其拟合曲线与原始曲线几乎重合，精度较高。2已知函数值346和边界条件：求三次样条插值函数并画出

15、其图形三次样条插值matlab程序如下：function ThrSample1 （x,y,y_1, y_N）syms t;if（length（x） = length（y） n = length（x）;x和y的维数不相等！end %维数检查A = diag（2*ones（1,n）; %求解m的系数矩阵u = zeros（n-2,1）;lamda = zeros（n-1,1）;c = zeros（n,1）; u（i-1） = （x（i）-x（i-1）/（x（i+1）-x（i-1）; lamda（i） = （x（i+1）-x（i）/（x（i+1）-x（i-1）; c（i） = 3*lamda（i）*

16、（y（i）-y（i-1）/（x（i）-x（i-1）+ . 3*u（i-1）*（y（i+1）-y（i）/（x（i+1）-x（i）; A（i, i+1） = u（i-1）; A（i, i-1） = lamda（i）; %形成系数矩阵及向量cc（1） = 2*y_1;c（n） = 2*y_N;m = followup（A,c）; %用追赶法求解方程组j=1;while j=n-1 h = x（j+1） - x（j）; f = y（j）*（2*（t-x（j）+h）*（t-x（j+1）2/h3 + . y（j+1）*（2*（x（j+1）-t）+h）*（t-x（j）2/h3 + . m（j）*（t-x（j

17、）*（x（j+1）-t）2/h2 - . m（j+1）*（x（j+1）-t）*（t-x（j）2/h2; fprintf（当x属于%d,%d,f=n, x（j）, x（j+1）; disp（f）; X=x（j）:x（j+1）; F=subs（f,t,X）; plot（X,F, hold on; j=j+1;function x=followup（A,b）n = rank（A）;for（i=1:Error: 对角有元素为0！end;d = ones（n,1）;a = ones（n-1,1）;c = ones（n-1）;n-1） a（i,1）=A（i+1,i）; c（i,1）=A（i,i+1）; d

18、（i,1）=A（i,i）;d（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）;x（n,1） = b（n,1）/d（n,1）;for（i=（n-1）:1） x（i,1） = （b（i,1）-c（i,1）*x（i+1,1）/d（i,1）;x=3 4 6;y=6 0 2;y_1=1;y_N=-1;ThrSample1 （x,y,y_1, y_N）由图像可以看出三次样条插值曲线的优越性，即光滑性好，收敛得到保证，局部性好。3观察物体的直线运动，得到如下数据时刻t0.91.93.03.95.0位移s10305180111求运动学方程。matlab程序如下：t=0 0.9 1.9 3.0 3.9 5.0;s=0 10 30 51 80 111;n=length（t）;a2=0; a2=t（i）2+a2;a3=0; a3=t（i）3+a3;a4=0; a4=t（i）4+a4;b1=0; b1=s（i）.*t（i）+b1;b2=0; b2=s（i）.*（t（i）.2）+b2;A=a2 a3;a3 a4b=b1;b2c=inv（A）*bx=0:0.01:5.0;y=c（1）.*x+c（