数值分析MATLAB编程.docx

题目1：用Gauss消去法解方程组x1+2x2+3x3=143x1+x2+5x3=202x1+5x2+2x3=18解：用MATLAB编写的程序如下：A=1 2 3;3 1 5;2 5 2;b=14 20 18;m,n=size(A);if rank(A)=rank(A,b)errorreturn;endc=n+1;A(:,c)=b;for k=1:n-1A(k+1:n,k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/A(k,k)*A(k,k:c);endx=zeros(length(b),1);x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);enddisp(x=);x=disp(x);题目2：已知4=2，9=3，用线性插值法求7的近似值。 x0=4:9; y0=x0.0.5; x=4:1:9; y=x.0.5; y2=interp1(x0,y0,x); figure(2),plot(x,y,b,x,y2,k:),grid disp(y2=)y2= disp(y2)2.0000 2.2361 2.4495 2.6458 2.8284 3.0000 题目三：设fx=ex，在-1,1上用Legendre多项式作fx的3次最佳平方逼近多项式。 a=-1;b=1; n=3; syms x; fx=exp(x); exp(x); for k=3F=xk*fx;d(k+1)=int(F,a,b);end d=reshape(d,n+1,1); H=hilb(n+1); a=Hda = -3696*exp(-1)+536*exp(1) 42900*exp(-1)-6180*exp(1) -105240*exp(-1)+15120*exp(1) 69300*exp(-1)-9940*exp(1) vpa(a)ans = 97.3166454843974255360582301282 -1016.953673622023958078606833033 2384.788857418174033934624995061 -1525.676101701956752911460193057题目四：用Gauss-Legendre求积公式（n=1,2）计算积分I=01x2exdx。 f=inline(1+t)/2)2*exp(1+t)/2)f = Inline function: f(t) = (1+t)/2)2*exp(1+t)/2) I1=1*feval(f,-0.57735)+1*feval(f,0.57735); I1=I1/2I1 =0.7119%n=1,两点Gauss-Legendre公式求解% f=inline(1+t)/2)2*exp(1+t)/2)f = Inline function: f(t) = (1+t)/2)2*exp(1+t)/2) x=0.7745966692,- 0.7745966692,0x = 0.7746 -0.7746 0 A=0.555555556, 0.555555556,0.888888889A = 0.5556 0.5556 0.8889 I2=0; for i=1:length(x)I2=I2+feval(f,x(i)/2*A(i);end I2I2 = 0.7183%n=2,三点Gauss-Legendre公式求解%题目五：用J法和GS法分别求解方程组10312-1031310x1x2x3=14-514其准确解为x*=1,1,1T。解：J法： A=10 3 1 14;2 -10 3 -5;1 3 10 14 MAXTIME=50; eps=1e-5; n,m=size(A); x=zeros(n,1); y=zeros(n,1); k=0; disp(X=); while 1 disp(x); for i=1:1:n s=0.0; for j=1:1:n if j=i s=s+A(i,j)*x(j); end y(i)=(A(i,n+1)-s)/A(i,i); end end for i=1:1:n maxeps=max(0,abs(x(i)-y(i); end if maxepsMAXTIME error return; end endGS法： format long;A=10 3 1 14;2 -10 3 -5;1 3 10 14 n,m=size(A); Maxtime=50; Eps=10E-5; x=zeros(1,n); disp(x=);x= for k=1:Maxtime disp(x); for i=1:n s=0.0; for j=1:n if i=j s=s+A(i,j)*x(j); end end x(i)=(A(i,n+1)-s)/A(i,i); end if sum(x-floor(x).2) f=inline(x+1)(1/3) disp(x=);x= x=feval(f,1.5); disp(x); 1.3572 i=1; while 1i=i+1; x=feval(f,x);disp(x);if i11;break;if x1E10;break;end;end;end;程序2： f=inline(x3-1) disp(x=);x= x=feval(f,1.5); disp(x); 2.3750 i=1; while 1i=i+1; x=feval(f,x);disp(x);if i11;break;if x1E10;break;end;end;end; 题目七：用幂法求矩阵A=110.5110.250.50.252的主特征值和主特征向量。 A=1 1 0.5;1 1 0.25;0.5 0.25 2A = 1.0000 1.0000 0.5000 1.0000 1.0000 0.2500 0.5000 0.2500 2.0000 Maxtime=100; Eps=1E-5; n=length(A); V=ones(n,1); k=0; m0=0; while k=Maxtime v=A*V; vmax,i=max(abs(v); m=v(i); V=v/m; if abs(m-m0) D=m; disp(D) 2.5365 disp(V) 0.7482 0.6497 1.0000题目八：分别取h=1,2,4，用经典R-K方法计算y=-y+x-e-1,y1=0,其准确解为yx=e-x+x-1-e-1。s=dsolve(Dy=-y+x-exp(-1),y(1)=0,x)clear;close;x=1:1:13;y=x-1-exp(-1)+exp(-x);disp(y) F=-y+x-exp(-1); a=1; b=13; h=1;%可通过改变h的值来满足题目要求% n=(b-a)/h; X=a:h:b; Y=zeros(1,n+1); Y(1)=0; for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=