矩阵分析与数值分析实验报告.docx
- 文档编号:9732277
- 上传时间:2023-02-06
- 格式:DOCX
- 页数:23
- 大小:566.44KB
矩阵分析与数值分析实验报告.docx
《矩阵分析与数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《矩阵分析与数值分析实验报告.docx(23页珍藏版)》请在冰豆网上搜索。
矩阵分析与数值分析实验报告
《矩阵分析与数值分析》实验报告
院系:
姓名:
学号:
所在班号:
任课老师:
一.设
,分别编制从小到大和从大到小的顺序程序计算
并指出有效位数。
程序如下:
functionsum3
j=input('请输入求和个数"j":
');
A=0;
B=0;
doubleB;
doubleA;
forn=2:
j
m=n^2-1;
t=1./m;
A=A+t;
end
disp('从小到大:
')
s=A
forn=j:
-1:
2
m=n^2-1;
t=1./m;
B=B+t;
end
disp('从大到小:
')
s=B
运行结果:
>>sum3
请输入求和个数"j":
100
从小到大:
s=0.740049504950495
从大到小:
s=0.740049504950495
>>sum3
请输入求和个数"j":
10000
从小到大:
s=0.749900004999506
从大到小:
s=0.749900004999500
>>sum3
请输入求和个数"j":
1000000
从小到大:
s=0.749999000000522
从大到小:
s=0.749999000000500
二、解线性方程组
1.分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组。
迭代法计算停止的条件为:
。
解:
(1)Jacobi迭代法程序代码:
functionjacobi(A,b,N)
clc;clear;
A=[-2100;1-210;01-21;001-2];
b=[-1000]';
N=100;
n=size(A,1);
D=diag(diag(A));
L=tril(-A,-1);
U=triu(-A,1);
Tj=inv(D)*(L+U);
cj=inv(D)*b;
tol=1e-06;
k=1;
formatlong
x=zeros(n,1);
whilek<=N
x(:
k+1)=Tj*x(:
k)+cj;
disp(k);disp('x=');disp(x(:
k+1));
ifnorm(x(:
k+1)-x(:
k)) disp('Theprocedurewassuccessful') disp('Condition||x^(k+1)-x^(k)|| disp(k);disp('x=');disp(x(: k+1)); break end k=k+1; end 结果输出 Theprocedurewassuccessful Condition||x^(k+1)-x^(k)|| 60 x= 0.79999879906731 0.59999842795870 0.39999805685009 0.199********505 (2)Gauss-Seidel迭代法程序代码: functiongauss_seidel(A,b,N) clc;clear; A=[-2100;1-210;01-21;001-2]; b=[-1000]'; N=100; n=size(A,1); D=diag(diag(A)); L=tril(-A,-1); U=triu(-A,1); Tg=inv(D-L)*U; cg=inv(D-L)*b; tol=1e-06; k=1; x=zeros(n,1); whilek<=N x(: k+1)=Tg*x(: k)+cg; disp(k);disp('x=');disp(x(: k+1)); ifnorm(x(: k+1)-x(: k)) disp('Theprocedurewassuccessful') disp('Condition||x^(k+1)-x^(k)|| disp(k);disp('x=');disp(x(: k+1)); break end k=k+1; end 结果输出 Theprocedurewassuccessful Condition||x^(k+1)-x^(k)|| 31 x= 0.79999921397935 0.59999897108561 0.39999916759077 0.199********539 2.用Gauss列主元消去法、QR方法求解如下方程组: (1)Gauss列主元消去法 程序代码: functionx=Gaussmain(A,b) clc;clear; formatlong A=[1212;253-2;-2-235;1323]; b=[47-10]'; N=length(A); x=zeros(N,1); y=zeros(N,1); c=0; d=0; A(: N+1)=b; fork=1: N-1 fori=k: 4 ifc d=i;c=abs(A(i,k)); end end y=A(k,: ); A(k,: )=A(d,: ); A(d,: )=y; fori=k+1: N c=A(i,k); forj=1: N+1 A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end end end b=A(: N+1); x(N)=b(N)/A(N,N); fork=N-1: -1: 1 x(k)=(b(k)-A(k,k+1: N)*x(k+1: N))/A(k,k); end 结果输出 ans= 18.00000000000000 -9.57142857142857 6.00000000000000 -0.42857142857143 (2)QR方法: 程序代码 functionQR(A,b) clc;clear; formatlong A=[1212;253-2;-2-235;1323]; b=[47-10]'; [Q,R]=qr(A) y=Q'*b; x=R\y 结果输出 Q= -0.31622776601684-0.04116934847963-0.751646028002830.57735026918962 -0.63245553203368-0.49403218175557-0.150********056-0.57735026918963 0.63245553203368-0.74104827263336-0.22549380840085-0.00000000000000 -0.31622776601684-0.452862833275940.601316822402260.57735026918963 R= -3.16227766016838-6.00832755431992-0.948683298050512.84604989415154 0-2.42899156029822-4.65213637819829-4.158********272 00-0.67648142520255-0.52615221960200 .0414******** x= 17.99999999999989 -9.57142857142851 5.99999999999997 -0.42857142857143 三、非线性方程的迭代解法 1.用Newton迭代法求方程 的根,计算停止的条件为: ; 编程如下: functionnewton(f,df,x,a,a0) symsx f=input('pleaseenteryourequation: ') a0=input('pleaseenteryoux(0): '); df=diff(f) e=1e-6; a1=a0+1; N=0; whileabs(a1-a0)>e a=a0-subs(f,a0)/subs(df,a0); a1=a0; a0=a; N=N+1; end fprintf('a=%0.6f',a) N 运行结果: >>newton pleaseenteryourequation: exp(x)+2^(-x)+2*cos(x)-6 f=exp(x)+2^(-x)+2*cos(x)-6 pleaseenteryoux(0): 2 df=exp(x)-2^(-x)*log (2)-2*sin(x) a=1.829384 N=4 2.利用Newton迭代法求多项式 的所有实零点,注意重根的问题。 由于不知道重根的个数,采用试探法求重根。 function[X]=newton2() clc;clear; symsx; delta=1e-06; f=inline('x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951'); df=inline('4*x^3-16.2*x^2+21.12*x-8.954'); u=inline('(x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951)/(4*x^3-16.2*x^2+21.12*x-8.954)'); du=inline('1-(x^4-27/5*x^3+264/25*x^2-4477/500*x+27951/10000)/(4*x^3-81/5*x^2+528/25*x-4477/500)^2*(12*x^2-162/5*x+528/25)'); j=1; fori=0: 1: 3 x0=i; while1 x1=x0-feval(u,x0)/feval(du,x0); err=abs(x1-x0); x0=x1; y=feval(f,x0); if(err X(j)=x1; j=j+1; break end end end X %运行结果: %X=1.10001.10002.10001.1000 四、数值积分 分别用复化梯形公式和Romberg公式计算积分 的近似值,要求误差不超过 ,并给出节点个数。 1)复化梯形公式 程序代码: function[I,step]=Trapezoid(f,a,b,eps) f='1/x';a=2;b=8;eps=1.0e-5; n=1; h=(b-a)/2; I1=2; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; whileabs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; fori=0: n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+... subs(sym(f),findsym(sym(f)),x1)); end end I=I2 step=n 结果输出 I= 1.38654458753674 step= 53 2)Romberg公式 程序代码: functionapprox=romberg(f,a,b,TOL) clc;clear; f=inline('1/x'); a=2;b=8; TOL=1e-05; A=zeros(20,20); A(1,1)=(b-a)*(feval(f,a)+feval(f,b))/2; h=(b-a)/2; A(2,1)=A(1,1)/2+h*feval(f,a+h); A(2,2)=(4*A(2,1)-A(1,1))/3; errest=abs(A(2,2)-A(1,1)/2); i=2; while(errest>TOL) i=i+1; h=h/2; sum=0.0; forj=1: 2: 2^(i-1)-1 sum=sum+feval(f,a+j*h); end; A(i,1)=A(i-1,1)/2+h*sum; forj=2: i power=4^(j-1); A(i,j)=(power*A(i,j-1)-A(i-1,j-1))/(power-1); end; errest=abs(A(i,i)-A(i-1,i-1))/2^(i-1); end; if(nargout==0) s=sprintf('\t\tapproximatevalueofintegral: \t%.12f\n',A(i,i)); s=sprintf('%s\t\terrorestimate: \t\t\t\t\t%.4e\n',s,errest); s=sprintf('%s\t\tnumberoffunctionevaluations: \t%d\n',s,2^(i-1)+1); disp(s) else approx=A(i,i); end 运行结果: approximatevalueofintegral: 1.386297441871 errorestimate: 8.7527e-006 numberoffunctionevaluations: 17 五、插值与逼近 1.给定 上的函数 ,请做如下的插值逼近: 构造等距节点分别取 , , 的Lagrange插值多项式; 构造分段线性取 的Lagrange插值多项式; 取Chebyshev多项式 的零点: , 作插值节点构造 的插值多项式 和上述的插值多项式均要求画出曲线图形(用不同的线型或颜色表示不同的曲线)。 程序代码 functionLagrange clc;clear;closeall; fori=1: 3; ifi==1 N=5; elseifi==2 N=8; else N=10; end f=inline('1/(1+25*x^2)'); x1=zeros(1,N+1); a=-1; b=1; fori=1: N+1 x1(i)=a+(i-1)*(b-a)/N; y1(i)=feval(f,x1(i)); end symsx ff=0; fori=1: N+1 f=1; forj=1: i-1 f=f*(x-x1(j))/(x1(i)-x1(j)); end forj=i+1: N+1 f=f*(x-x1(j))/(x1(i)-x1(j)); end ff=f*y1(i)+ff; f=1; end ff=collect(ff,x); ff=vpa(ff,4); y=ff; p=ezplot(y,[a,b]);grid YLIM([-0.10.6]); ifN==5 set(p,'Color','black'); set(p,'LineStyle','--'); lagrange_5=y elseifN==8 set(p,'Color','r'); set(p,'LineStyle','--'); lagrange_8=y else set(p,'Color','b'); set(p,'LineStyle','--') lagrange_10=y end holdon; xlabel('x');ylabel('y'); title('y=p(x)');holdon; end Lag_Cheb(); x=-1: 0.01: 1; y=1./(1+25*x.^2); acu=plot(x,y);gridon;holdon set(acu,'Color','m'); set(acu,'LineStyle','-'); legend('N=5','N=8','N=10','Cheb,N=10','¾«È·Öµ'); functionLag_Cheb() f=inline('1/(1+25*x^2)'); N=10; x1=zeros(1,N+1); a=-1; b=1; fori=1: N+1 x1(i)=cos((2*i-1)*pi/(2*N)); y1(i)=feval(f,x1(i)); end symsx ff=0; fori=1: N+1 f=1; forj=1: i-1 f=f*(x-x1(j))/(x1(i)-x1(j)); end forj=i+1: N+1 f=f*(x-x1(j))/(x1(i)-x1(j)); end ff=f*y1(i)+ff; f=1; end ff=collect(ff,x); ff=vpa(ff,4); yy=ff; ff=collect(ff,x); yy=ff; lagrange_chebshev_10=yy cheb=ezplot(yy,[a,b]);gridon YLIM([-0.10.6]); set(cheb,'Color','g'); set(cheb,'LineStyle',': '); 结果输出 lagrange_5=.5673+1.202*x^4-1.731*x^2 lagrange_8=1.+53.69*x^8-102.8*x^6+61.37*x^4-13.20*x^2 lagrange_10=1.-220.9*x^10+494.9*x^8-381.4*x^6+123.4*x^4-16.86*x^2 lagrange_chebshev_10=.7413-5.359*x^10-.4e-2*x^9+18.96*x^8-.1321*x^7- 25.78*x^6+.10*x^5+16.81*x^4+.5e-2*x^3-5.336*x^2+.5288e-3*x 2.已知函数值 0 1 2 3 4 5 6 7 8 9 10 2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80 和边界条件: . 求三次样条插值函数 并画出其图形. 程序代码: functionSpline3_1(x,y,df0,dfn) formatshort x=[012345678910]; y=[2.513.304.044.705.225.545.785.405.575.705.80]; plot(x,y,'g*','MarkerSize',15);holdon; df0=1; dfn=0; n=length(x); h=zeros(n-1,1); lan=zeros(n-2,1); mu=zeros(n-2,1); g=zeros(n-2,1); m=zeros(n,1); m (1)=df0; m(n)=dfn; fori=1: n-1 h(i)=x(i+1)-x(i); end fori=1: n-2 mu(i)=h(i)/(h(i+1)+h(i)); lan(i)=h(i+1)/(h(i+1)+h(i)); g(i)=3*(mu(i)*(y(i+2)-y(i+1))/h(i+1)+lan(i)*(y(i+1)-y(i))/h(i)); end A=zeros(n-2,n-2); A(1,1)=2; A(1,2)=mu (1); A(n-2,n-2)=lan(n-2); fori=2: n-2 A(i,i)=2; A(i,i-1)=mu(i); A(i-1,i)=lan(i); end g (1)=g (1)-lan (1)*df0; g(n-2)=g(n-2)-mu(n-2)*dfn; b=A\g; fori=2: n-1 m(i)=b(i-1); end symsz fori=1: n-1 xx=x(i): 0.01: x(i+1); sx1=y(i)*(xx-x(i+1)).^2.*(h(i)+2*(xx-x(i)))/h(i)^3; sx2=y(i+1)*(xx-x(i)).^2.*(h(i)-2*(xx-x(i+1)))/h(i)^3; sx3=m(i)*(xx-x(i+1)).^2.*(xx-x(i))/h(i)^2; sx4=m(i+1)*(xx-x(i)).^2.*(xx-x(i+1))/h(i)^2; sx=sx1+sx2+sx3+sx4; z1=y(i)*(z-x(i+1)).^2.*(h(i)+2*(z-x(i)))/h(i)^3; z2=y(i+1)*(z-x(i)).^2.*(h(i)-2*(z-x(i+1)))/h(i)^3; z3=m(i)*(z-x(i+1)).^2.*(z-x(i))/h(i)^2; z4=m(i+1)*(z-x(i)).^2.*(z-x(i+1))/h(i)^2; zz=z1+z2+z3+z4; zz=collect(zz); vpa(zz,4) plot(xx,sx,'r'); holdon; gridon; end 程序输出 ans=2.510+.1379*z^3-.3479*z^2+z ans=2.692-.4369e-1*z^3+.1969*z^2+.4552*z ans=2.287+.6872e-2*z^3-.1065*z^2+1.062*z ans=3.655-.4380e-1*z^3+.3495*z^2-.3060*z ans=-6.080+.1083*z^3-1.476*z^2+6.995*z ans=41.14-.2694*z^3+4.190*z^2-21.34*z ans=-109.8+.4295*z^3-8.390*z^2+54.14*z ans=133.0-.2784*z^3+6.475*z^2-49.91*z ans=-57.75+.9411e-1*
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 分析 数值 实验 报告