同濟大學數(shù)值分析matlab編程題匯編

1、MATLAB 編程題庫1.下面的數(shù)據(jù)表近似地滿足函數(shù)，請適當變換成為線性最小二乘問題，編程求最好的系數(shù)，并在同一個圖上畫出所有數(shù)據(jù)和函數(shù)圖像.解：>> x=-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995'>> y=0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625'>> A=x ones(8,1 -x.2.*y;>> z=Ay;>> a=z(1; b=z(2; c=z(3;>> xh=-1:0.1:1'

2、;>> yh=(a.*xh+b./(1+c.*xh.2;>> plot(x,y,'r+',xh,yh,'b*'2.若在Matlab工作目錄下已經(jīng)有如下兩個函數(shù)文件，寫一個割線法程序，求出這兩個函數(shù)精度為的近似根，并寫出調(diào)用方式：文件一文件二function v = f(xv = x .* log(x - 1;function z = g(yz = y.5 + y - 1;解：>> edit gexianfa.mfunction x iter=gexianfa(f,x0,x1,toliter=0;x=x1;while(abs(f

3、eval(f,x>toliter=iter+1;x=x1-feval(f,x1.*(x1-x0./(feval(f,x1-feval(f,x0;x0=x1;x1=x;end>> edit f.mfunction v=f(xv=x.*log(x-1;>> edit g.mfunction z=g(yz=y.5+y-1;>> x1 iter1=gexianfa('f',1,3,1e-10x1 =1.7632iter1 =6>> x2 iter2=gexianfa('g',0,1,1e-10x2 =0.7549it

4、er2 =83.使用GS迭代求解下述線性代數(shù)方程組：解：>> edit gsdiedai.mfunction x iter=gsdiedai(A,x0,b,tolD=diag(diag(A;L=D-tril(A;U=D-triu(A;iter=0;x=x0;while(norm(b-A*x./norm(b>toliter=iter+1;x0=x;x=(D-L(U*x0+b;end>> A=5 2 1;-1 4 2;1 -3 10;>> b=-12 10 3'>> tol=1e-4;>> x0=0 0 0'>

5、> x iter=gsdiedai(A,x0,b,tol;>> xx =-3.09101.23720.9802>> iteriter =64.用四階Range-kutta方法求解下述常微分方程初值問題（取步長h=0.01）解：>> edit ksf2.mfunction v=ksf2(x,yv=y+exp(x+x.*y;>> a=1;b=2;h=0.01;>> n=(b-a./h;>> x=1:0.01:2;>> y(1=2;>> for i=2:(n+1k1=h*ksf2(x(i-1,y(i

6、-1;k2=h*ksf2(x(i-1+0.5*h,y(i-1+0.5*k1;k3=h*ksf2(x(i-1+0.5*h,y(i-1+0.5*k2;k4=h*ksf2(x(i-1+h,y(i-1+k3;y(i=y(i-1+(k1+2*k2+2*k3+k4./6;end>> y調(diào)用函數(shù)方法>> edit Rangekutta.mfunction x y=Rangekutta(f,a,b,h,y0x=a:h:b;n=(b-a/h;y(1=y0;for i=2:(n+1k1=h*(feval(f,x(i-1,y(i-1;k2=h*(feval(f,x(i-1+0.5*h,y(i

7、-1+0.5*k1;k3=h*(feval(f,x(i-1+0.5*h,y(i-1+0.5*k2;k4=h*(feval(f,x(i-1+h,y(i-1+k3;y(i=y(i-1+(k1+2*k2+2*k3+k4./6;end>> x y=Rangekutta('ksf2',1,2,0.01,2;>> y5.取，請編寫Matlab程序，分別用歐拉方法、改進歐拉方法在上求解初值問題。解：>> edit Euler.mfunction x y=Euler(f,a,b,h,y0x=a:h:b;n=(b-a./h;y(1=y0;for i=2:(n+

8、1y(i=y(i-1+h*feval(f,x(i-1,y(i-1;end>> edit gaijinEuler.mfunctionx y=gaijinEuler(f,a,b,h,y0x=a:h:b;n=(b-a./h;y(1=y0;for i=2:(n+1y1=y(i-1+h*feval(f,x(i-1,y(i-1;y2=y(i-1+h*feval(f,x(i,y1;y(i=(y1+y2./2;end>> edit ksf3.mfunction v=ksf3(x,yv=x.3-y./x;>> x y=Euler('ksf3',1,2,0.2,

9、0.4x =1.0000 1.2000 1.4000 1.6000 1.8000 2.0000y =0.4000 0.5200 0.7789 1.2165 1.8836 2.8407>> x y=gaijinEuler('ksf3',1,2,0.2,0.4x =1.0000 1.2000 1.4000 1.6000 1.8000 2.0000y =0.4000 0.5895 0.9278 1.4615 2.2464 3.34666.請編寫復合梯形積分公式的Matlab程序，計算下面積分的近似值，區(qū)間等分。編寫辛普森積分公式的Matlab程序，計算下面積分的近似值，區(qū)

10、間等分。、解：>> edit tixingjifen.mfunction s=tixingjifen(f,a,b,nx=linspace(a,b,(n+1;y=zeros(1,length(x;y=feval(f,xh=(b-a./n;s=0.5*h*(y(1+2*sum(y(2:n+y(n+1;end>> edit simpson.mfunction I=simpson(f,a,b,nh=(b-a/n;x=linspace(a,b,2*n+1;y=feval(f,x;I=(h/6*(y(1+2*sum(y(3:2:2*n-1+4*sum(y(2:2:2*n+y(2*n

11、+1;>> edit ksf4.mfunction v=ksf4(xv=1./(x.2+1;>> tixingjifen('ksf4',0,1,20ans =0.7853>>simpson('ksf4',0,1,10ans =0.7854>> edit ksf5.mfunction v=ksf5(xif(x=0v=1;elsev=sin(x./x;end(第二個函數(shù)ksf5調(diào)用求積函數(shù)時，總顯示有錯誤：“NaN”，還沒調(diào)試好。見諒！7.用迭代方法對下面方程組求解，取初始向量。解：>>edit Jacob

12、i.mfunctionx iter=Jacobi(A,x0,b,tolD=diag(diag(A;L=D-tril(A;U=D-triu(A;x=x0;iter=0;while(norm(A*x-b/norm(b>toliter=iter+1;x0=x;x=D(L+U*x0+b;end>> A=2 4 -4;3 3 3;4 4 2;>> b=2 -3 -2'>> x0=3 2 -1'>> x,iter=Jacobi(A,x0,b,1e-4x =1-1-1iter =38.用牛頓法求解方程在附近的根。解：>> ed

13、it Newton.mfunction x iter=Newton(f,g,x0,toliter=0;x=x0;while(abs(feval(f,x>tolx0=x;x=x0-feval(f,x0./feval(g,x0;iter=iter+1;end>> edit ksf6.mfunction v=ksf6(xv=x*cos(x+2;>> edit ksg6.mfunction z=ksg(yz=y.5+y-1;>> x iter=Newton('ksf6','ksg6',2,1e-4x =2.4988iter =3

14、9.分別用改進乘冪法、反冪法計算矩陣A的按模最大特征值及其對應的特征向量、按模最小特征值及其對應的特征向量。解：>> edit ep.mfunction t,x=ep(A,x0,toltv0 ti0=max(abs(x0;lam0=x0(ti0;x0=x0./lam0;x1=A*x0;tv1 ti1=max(abs(x1;lam1=x1(ti1;x1=x1./lam1;while(abs(lam0-lam1>tolx0=x1;lam0=lam1;x1=A*x0;tv1 ti1=max(abs(x1;lam1=x1(ti1;x1=x1./lam1;endt=lam1;x=x1;

15、>> edit fanep.mfunctiont,x=fanep(A,x0,toltv0 ti0=max(abs(x0;lam0=x0(ti0;x0=x0./lam0;x1=Ax0;tv1 ti1=max(abs(x1;lam1=x1(ti1;x1=x1./lam1;while(abs(1/lam0-1/lam1>tolx0=x1;lam0=lam1;x1=Ax0;tv1 ti1=max(abs(x1;lam1=x1(ti1;x1=x1./lam1;endt=1/lam1;x=x1;>> A=12 6 -6;6 16 2;-6 2 16;>>x0=1

16、0.5 -0.5'>>tol=1e-4;>>t,x=ep(A,x0,tolt =21.5440x =1.00000.7953-0.7953>> A=12 6 -6;6 16 2;-6 2 16;>>x0=1 0.5 -0.5'>>tol=1e-4;t,x=fanep(A,x0,tolt =4.4560x =1.0000-0.62870.628710.將積分區(qū)間n等分，用復合梯形求積公式計算定積分,比較不同值時的誤差（畫出平面上 log(n-log(Error圖）解：>> edit ksf7.mfunctio

17、n v=ksf7(xv=sqrt(1+x.2;>> I=quad('ksf7',1,3;>> n=1:100;>> for i=1:100x(i=tixingjifen('ksf7',1,3,i;error(i=abs(I-x(i;end>> plot(log10(n,log10(error(n11. 用迭代方法對下面方程組求解，取初始向量。解：>> edit sor.mfunction x,iter=sor(A,x0,b,omega,tolD=diag(diag(A;L=D-tril(A;U=D-tr

18、iu(A;x=x0;iter=0;while(norm(A*x-b/norm(b>toliter=iter+1;x0=x;x=(D-omega*L(omega*b+(1-omega*D*x+omega*U*x;end>> A=2 -1 0;-1 3 -1;0 -1 2;>> b=1 8 -5'>> x0=1 0 -1'>> x,iter=sor(A,x0,b,1.1,1e-4x =1.99993.0000-1.0000iter =5其他>> A=1 2 3 4;5 6 7 8;>> AA =1 2 3 45 6 7 8>> m n=size(Am =2n =4>> x=1 2 3 4 5;>> xx =1 2 3 4 5>> length(xans =5>> A=2 3 4;1 1 9; 1 2 -6;>> AA =2 3 41 1 91 2 -6>> L U=lu(A;>> LL =1.0000 0 00.5000 1.0000 00.5000 -1.0000 1.0000>> UU =2.0000 3.0000 4.000