同济大学数值分析matlab编程题汇编

MATLAB编程题库1.下面的数据表近似地满足函数，请适当变换成为线性最小二乘问题，编程求最好的系数，并在同一个图上画出所有数据和函数图像.解：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工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数精度为的近似根，并写出调用方式：文件一文件二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程序，分别用欧拉方法、改进欧拉方法在上求解初值问题。解：>>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.请编写复合梯形积分公式的Matlab程序，计算下面积分的近似值，区间等分。编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分。、解：>>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(第二个函数‘ksf5’调用求积函数时，总显示有错误：“NaN”，还没调试好。见谅！)

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.分别用改进乘幂法、反幂法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。解：>>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.将积分区间n等分，用复合梯形求积公式计算定积分,比较不同值时的误差（画出平面上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.编程实现求解满足下列条件的区间[-1,2]上的三次样条函数S(x),并画出此样条函数的图形: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.编程实现求解满足下列条件的区间[-1,2]上的三次样条函数S(x),并画出此样条函数的图形:xi-1012f(xi)-1010f(xi)'0-1functionsplx=[-1012