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

MATLAB編程題庫1.下面的數(shù)據(jù)表近似地滿足函數(shù)，請適當(dāng)變換成為線性最小二乘問題，編程求最好的系數(shù)，并在同一個圖上畫出所有數(shù)據(jù)和函數(shù)圖像.解：x=[-0.931-0.586-0.362-0.2130.0080.5440.6280.995]';y=[0.3560.6060.6870.8020.8230.8010.7180.625]';A=[xones(8,1)-x.^2.*y];z=A\y;a=z(1);b=z(2);c=z(3);xh=-1:0.1:1;yh=(a.*xh+b)./(1+c.*xh.^2);plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目錄下已經(jīng)有如下兩個函數(shù)文件，寫一個割線法程序，求出這兩個函數(shù)精度為的近似根，并寫出調(diào)用方式：文件一文件二functionv=f(x)v=x.*log(x)-1;functionz=g(y)z=y.^5+y-1;解：>>editgexianfa.mfunction[xiter]=gexianfa(f,x0,x1,tol)iter=0;while(norm(x1-x0)>tol)iter=iter+1;x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0));x0=x1;x1=x;end>>editf.mfunctionv=f(x)v=x.*log(x)-1;>>editg.mfunctionz=g(y)z=y.^5+y-1;>>[x1iter1]=gexianfa('f',1,3,1e-10)x1=1.7632iter1=6>>[x2iter2]=gexianfa('g',0,1,1e-10)x2=0.7549iter2=8

5.取，請編寫Matlab程序，分別用歐拉方法、改進(jìn)歐拉方法在上求解初值問題。解：>>editEuler.mfunction[xy]=Euler(f,a,b,h,y0)x=[a:h:b];n=(b-a)./h;y(1)=y0;fori=2:(n+1)y(i)=y(i-1)+h*feval(f,x(i-1),y(i-1));end>>editgaijinEuler.mfunction[xy]=gaijinEuler(f,a,b,h,y0)x=[a:h:b];n=(b-a)./h;y(1)=y0;fori=2:(n+1)y1=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>>editksf3.mfunctionv=ksf3(x,y)v=x.^3-y./x;>>[xy]=Euler('ksf3',1,2,0.2,0.4)x=1.00001.20001.40001.60001.80002.0000y=0.40000.52000.77891.21651.88362.8407>>[xy]=gaijinEuler('ksf3',1,2,0.2,0.4)x=1.00001.20001.40001.60001.80002.0000y=0.40000.58950.92781.46152.24643.3466

6.請編寫復(fù)合梯形積分公式的Matlab程序，計算下面積分的近似值，區(qū)間等分。編寫辛普森積分公式的Matlab程序，計算下面積分的近似值，區(qū)間等分。、解：>>edittixingjifen.mfunctions=tixingjifen(f,a,b,n)x=linspace(a,b,(n+1));y=zeros(1,length(x));y=feval(f,x)h=(b-a)./n;s=0.5*h*(y(1)+2*sum(y(2:n))+y(n+1));end>>editsimpson.mfunctionI=simpson(f,a,b,n)h=(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+1));>>editksf4.mfunctionv=ksf4(x)v=1./(x.^2+1);>>tixingjifen('ksf4',0,1,20)ans=0.7853>>simpson('ksf4',0,1,10)ans=0.7854>>editksf5.mfunctionv=ksf5(x)if(x==0)v=1;elsev=sin(x)./x;end(第二個函數(shù)‘ksf5’調(diào)用求積函數(shù)時，總顯示有錯誤：“NaN”，還沒調(diào)試好。見諒！)

7.用迭代方法對下面方程組求解，取初始向量。解：>>editJacobi.mfunction[xiter]=Jacobi(A,x0,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=x0;iter=0;while(norm(A*x-b)/norm(b)>tol)iter=iter+1;x0=x;x=D\((L+U)*x0+b);end>>A=[24-4;333;442];>>b=[2-3-2]';>>x0=[32-1]';>>[x,iter]=Jacobi(A,x0,b,1e-4)x=1-1-1iter=3

8.用牛頓法求解方程在附近的根。解：>>editNewton.mfunction[xiter]=Newton(f,g,x0,tol)iter=0;done=0while~donex=x0-feval(f,x0)/feval(g,x0);done=norm(x-x0)<=tol;iter=iter+1;if~done,x0=x;endend>>editksf6.mfunctionv=ksf6(x)v=x*cos(x)+2;>>editksg6.mfunctionz=ksg(y)z=y.^5+y-1;>>[xiter]=Newton('ksf6','ksg6',2,1e-4)x=2.4988iter=3

9.分別用改進(jìn)乘冪法、反冪法計算矩陣A的按模最大特征值及其對應(yīng)的特征向量、按模最小特征值及其對應(yīng)的特征向量。解：>>editep.mfunction[t,x]=ep(A,x0,tol)[tv0ti0]=max(abs(x0));lam0=x0(ti0);x0=x0./lam0;x1=A*x0;[tv1ti1]=max(abs(x1));lam1=x1(ti1);x1=x1./lam1;while(abs(lam0-lam1)>tol)x0=x1;lam0=lam1;x1=A*x0;[tv1ti1]=max(abs(x1));lam1=x1(ti1);x1=x1./lam1;endt=lam1;x=x1;>>editfanep.mfunction[t,x]=fanep(A,x0,tol)[tv0ti0]=max(abs(x0));lam0=x0(ti0);x0=x0./lam0;x1=A\x0;[tv1ti1]=max(abs(x1));lam1=x1(ti1);x1=x1./lam1;while(abs(1/lam0-1/lam1)>tol)x0=x1;lam0=lam1;x1=A\x0;[tv1ti1]=max(abs(x1));lam1=x1(ti1);x1=x1./lam1;endt=1/lam1;x=x1;>>A=[126-6;6162;-6216];>>x0=[10.5-0.5]';>>tol=1e-4;>>[t,x]=ep(A,x0,tol)t=21.5440x=1.00000.7953-0.7953>>A=[126-6;6162;-6216];>>x0=[10.5-0.5]';>>tol=1e-4;[t,x]=fanep(A,x0,tol)t=4.4560x=1.0000-0.62870.6287

10.將積分區(qū)間n等分，用復(fù)合梯形求積公式計算定積分,比較不同值時的誤差（畫出平面上log(n)-log(Error)圖）．解：>>editksf7.mfunctionv=ksf7(x)v=sqrt(1+x.^2);>>I=quad('ksf7',1,3);>>n=1:100;>>fori=1:100x(i)=tixingjifen('ksf7',1,3,i);error(i)=abs(I-x(i));end>>plot(log10(n),log10(error(n)))

11.用迭代方法對下面方程組求解，取初始向量。解：>>editsor.mfunction[x,iter]=sor(A,x0,b,omega,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=x0;iter=0;while(norm(A*x-b)/norm(b)>tol)iter=iter+1;x0=x;x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x);end>>A=[2-10;-13-1;0-12];>>b=[18-5]';>>x0=[10-1]';>>[x,iter]=sor(A,x0,b,1.1,1e-4)x=1.99993.0000-1.0000iter=512.編程實現(xiàn)求解滿足下列條件的區(qū)間[-1,2]上的三次樣條函數(shù)S(x),并畫出此樣條函數(shù)的圖形:xi-1012f(xi)-1010f(xi)'0-1functionsplx=[-1012]y=[0-1010-1]pp=csape(x,y,'complete')[breaks,coefs,npolys,ncoefs,dim]=unmkpp(pp)xh=-1:0.1:2if-1<=xh<=0yh=coefs(1,1)*(xh+1).^3+coefs(1,2)*(xh+1).^2+coefs(1,3)*(xh+1)+coefs(1,4)elseif0<=xh<=1yh=coefs(2,1)*(xh).^3+coefs(2,2)*(xh).^2+coefs(2,3)*(xh)+coefs(2,4)elseif1<=xh<=2yh=coefs(3,1)*(xh-1).^3+coefs(3,2)*(xh-1).^2+coefs(3,3)*(xh-1)+coefs(3,4)elsereturnendplot(xh,yh,'r+')13.二分法程序functionx=bisect(f,a,b,tol)ifnargin<4tol=1e-12endfa=feval(f,a)fb=feval(f,b)whileabs(a-b)>tolx=(a+b)/2fx=feval(f,x)ifsign(fx)==sign(fa)a=xfa=fxelseifsign(fx)==sign(fb)b=xfb=fxelsereturnendend

其他>>A=[1234;5678];>>AA=12345678>>[mn]=size(A)m=2n=4>>x=[12345];>>xx=12345>>length(x)ans=5>>A=[234;119;12-6];>>AA=23411912-6>>[LU]=lu(A);>>LL=1.0000000.50001.000000.5000-1.00001.0000>>UU=2.00003.00004.00000-0.50007.000000-1.0000>>x=linspace(0,1,11);>>x'ans=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000>>x=[1234];>>y=[6111827];>>p=polyfit(x,y,2)p=1.00002.00003.0000>>diag(ones(4,1),1)ans=010000010000010000010000012.編程實現(xiàn)求解滿足下列條件的區(qū)間[-1,2]上的三次樣條函數(shù)S(x),并畫出此樣條函數(shù)的圖形:xi-1012f(xi)-1010f(xi)'0-1functionsplx=[-1012