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

1、MATLAB编程题库ax+bi下面的数据表近似地满足函数y二，请适当变换成为线性最小二乘问题，编程求1+cx2最好的系数a,b,c，并在同一个图上画出所有数据和函数图像.xi-0.931-0.586-0.362-0.2130.0080.5440.6280.995yi0.3560.6060.6870.8020.8230.8010.7180.625解：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.A2.*y;z=Ay;a

2、=z(1);b=z(2);c=z(3);xh=-1:0.1:1;yh=(a.*xh+b)./(1+c.*xh.A2);plot(x,y,'r+',xh,yh,'b*')182.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数精度为1010的近似根，并写出调用方式：文件一文件二functionv=f(x)functionz=g(y)v=x.*log(x)-1;z=y.A5+y-1;解：>>editgexianfa.mfunctionxiter=gexianfa(f,x0,x1,tol)iter=0;while(norm

3、(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=y45+y-l;>>x1iter1=gexianfa('f',1,3,1e-10)xl=l.7632iterl=6>>x2iter2=gexianfa('g',0,l,le-l0)x2=0.7549iter2=

4、83.使用GS迭代求解下述线性代数方程组：5x2xx12123-一x+4x+2x10123x3x+10x=31123解：>>editgsdiedai.mfunctionxiter=gsdiedai(A,x0,b,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);iter=0;x=x0;while(norm(b-A*x)./norm(b)>tol)iter=iter+1;x0=x;x=(D-L)(U*x0+b);end>>A=521;-142;1-310;>>b=-12103'>>tol=1e-4;

5、>>x0=000'>>xiter=gsdiedai(A,x0,b,tol);>>xx=-3.09101.23720.9802>>iteriter=64.用四阶Range-kutta方法求解下述常微分方程初值问题(取步长h=0.01)dy=y+ex+xy,-dxy(1)=2解：>>editksf2.mfunctionv=ksf2(x,y)v=y+exp(x)+x.*y;>>a=1;b=2;h=0.01;>>n=(b-a)./h;>>x=1:0.01:2;>>y(1)=2;>&

6、gt;fori=2:(n+1)k1=h*ksf2(x(i-1),y(i-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调用函数方法>>editRangekutta.mfunctionxy=Rangekutta(f,a,b,h,y0)x=a:h:b;n=(b-a)/h;y(1)=y0;fori=2:(n+1)k1=h*(fev

7、al(f,x(i-1),y(i-1);k2=h*(feval(f,x(i-1)+0.5*h,y(i-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>>xy=Rangekutta('ksf2',1,2,0.01,2);>>y5.取h=0.2，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在1<x<2上求解初值问题。y(1)=0.4解：>

8、>editEuler.mfunctionxy=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.mfunctionxy=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;en

9、d>>editksf3.mfunctionv=ksf3(x,y)v=x43-y./x;>>xy=Euler('ksf3',1,2,0.2,0.4)x=1.60001.80002.00001.00001.20001.4000y=0.40000.52000.77891.21651.88362.8407>>xy=gaijinEuler('ksf3',1,2,0.2,0.4)x=1.60001.80002.00001.00001.20001.4000y=0.40000.58950.92781.46152.24643.34666.请编写

10、复合梯形积分公式的Matlab程序，计算下面积分的近似值，区间等分n=20。编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分n=10。1X2+1dx、p1sinx,dxJ-1x解：>>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

11、,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./(x42+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

12、;elsev=sin(x)./x;end(第二个函数ksf5'调用求积函数时，总显示有错误：“NaN”，还没调试好。见谅！)7.用Jacobi迭代方法对下面方程组求解，取初始向量X(0)=(3,2,-1)T。243344-4_(X)r213X=-3221XJ3<-2丿解：>>editJacobi.mfunctionxiter=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(

13、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=8.用牛顿法求解方程xcosx+2=0在x0二2附近的根。解：>>editNewton.mfunctionxiter=Newton(f,g,x0,tol)iter=0;done=0whiledonex=x0-feval(f,x0)/feval(g,x0);done=norm(x-x0)<=tol;iter=iter+1;ifdone,

14、x0=x;endend>>editksf6.mfunctionv=ksf6(x)v=x*cos(x)+2;>>editksg6.mfunctionz=ksg(y)z=y45+y-l;>>xiter=Newton('ksf6','ksg6',2,1e-4)x=2.4988iter=39分别用改进乘幕法、反幕法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。126-6A=6162-6216解：>>editep.mfunctiont,x=ep(A,x0,tol)tv0ti0=max(abs

15、(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.mfunctiont,x=fanep(A,x0,tol)tv0ti0=max(abs(x0);lam0=x0(ti0);x0=x0./lam0;x1=Ax0;tv1ti1=max

16、(abs(x1);lam1=x1(ti1);x1=x1./lam1;while(abs(1/lam0-1/lam1)>tol)x0=x1;lam0=lam1;x1=Ax0;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=>>A=126-6;6162;-6216;>>x0=10.5-0.5'>&

17、gt;tol=1e-4;t,x=fanep(A,x0,tol)t=4.4560x=1.00000.7953-0.795321.5440x=1.0000-0.62870.628710.将积分区间n等分，用复合梯形求积公式计算定积分I'i+x2dx，比较不同n值时的误1差(画出平面上log(n)-log(Error)图).解：>>editksf7.mfunctionv=ksf7(x)v=sqrt(1+x.A2);>>I=quad('ksf7'，1，3);>>n=1:100;>>fori=1:100x(i)=tixingjife

18、n('ksf7'，1，3，i);error(i)=abs(I-x(i);end>>plot(log10(n)，log10(error(n)11.用SOR迭代方法对下面方程组求解，取初始向量X(0)二(1,0,-1)t。2-10-10_(X)1Xr13-1=82-121X丿3<-5丿解：>>editsor.mfunctionx,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)ite

19、r=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=-1012y=0-1010-1pp=cs

20、ape(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).A3+coefs(2,2)*(xh).A2+coefs(2,3)*(xh)+coefs(2,4)elseif1<=xh<=2yh=coefs(3,1)*(xh-1).A3+

21、coefs(3,2)*(xh-1).A2+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=

22、1234;5678;>>AA=123567>>mn=size(A)m=2n=4>>x=12345;>>xx=123>>length(x)ans=545>>A=234;119;12-6;>>AA=23411912-6>>LU=lu(A);>>LL=1.000000.50001.00000.5000-1.0000>>UU=2.00003.00000-0.500000001.00004.00007.0000-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&