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

1、MATLAB编程题库ax b 1.下面的数据表近似地满足函数y 2 ,请适当变换成为线性最小二乘问题,编程求1 cx最好的系数a,b,c,并在同一个图上画出所有数据和函数图像.xi0.9310.5860.3620.213 0.008 0.544 0.628 0.995yi0.3560.6060.6870.8020.823 0.801 0.718 0.625解：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 one

2、s(8,1) -x.A2.*y;z=Ay;a=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*')2.假设在Matlab工作目录下已经有如下两个函数文件,写一个割线法程序,求出这两个函数 精度为10 10的近似根,并写出调用方式：文件一文件二function v = f(x)v = x .* log(x) - 1;function z = g(y) z = y.A5 + y - 1;解：> > edit gexianfa.mfunct

3、ion x iter=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> > edit f.m function v=f(x) v=x.*log(x)-1;> > edit g.m function z=g(y) z=y.A5+y-1;> > x1 iter1=gexianfa('f,1,3,1e-10)x1 =1.7632iterl =6

4、>> x2 iter2=gexianfa('g',0,1,1e-10)x2 =0.7549iter2 =83.使用GS迭代求解下述线性代数方程组:15xi + 2x2 + X3 = - 12于 Xi + 4X2 + 2x3 = 10 羯 3x2+ 10x3 = 3解：> > edit gsdiedai.mfunction x iter=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=i

5、ter+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'> > x iter=gsdiedai(A,x0,b,tol);>>xx =-3.09101.23720.9802>>iteriter =6h=0.01)4用四阶Range-kutta方法求解下述常微分方程初值问题取步长2dyx,?= y+ e + xy,另dx?y(1)= 2解：>> edit

6、ksf2.m function v=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;>>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+k

7、4)./6;end>>y调用函数方法>> edit Rangekutta.mfunction x y=Rangekutta(f,a,b,h,y0)x=a:h:b;n=(b-a)/h;y(1)=y0;fori=2:(n+1)k1=h*(feval(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).

8、/6;end>> x y=Rangekutta('ksf2',1,2,0.01,2);>>y5.取h 0.2,请编写Matlab程序,分别用欧拉方法、改进欧拉方法在1 x问题.2上求解初值?dy=x3- y, idx x ?y(1)= 0.4解：>> edit Euler.mfunction x y=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>> edit gaijinEuler.mf

9、unctionx y=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).end>> edit ksf3.mfunction v=ksf3(x,y)v=x.A3-y./x;>>x y=Euler('ksf3',1,2,0.2,0.4)x =1.00001.20001.40001.60001.80002.0000y =0.40000.5

10、2000.77891.21651.88362.8407>>x y=gaijinEuler('ksf3',1,2,0.2,0.4)x =1.00001.20001.40001.60001.80002.0000y =0.40000.58950.92781.46152.24643.34666.请编写复合梯形积分公式的Matlab程序,计算下面积分的近似值,区间等分n 20.编写辛普森积分公式的Matlab程序,计算下面积分的近似值,区间等分n 10.1,dx、 0 皿 dx0 x2 + 1 Ol x解：> > edit tixingjifen.mfunctio

11、n s=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> > edit simpson.mfunction I=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);> > edit ksf4.

12、mfunction v=ksf4(x)v=1./(x.A2+1);>>tixingjifen('ksf4',0,1,20) ans =0.7853>>simpson('ksf4',0,1,10) ans =0.7854> > edit ksf5.mfunction v=ksf5(x)if(x=0)v=1;elsev=sin(x)./x;end(第二个函数ksf调用求积函数时,总显示有错误：“NaN；还没调试好.见谅！)7.用Jacobi迭代方法对下面方程组求解,取初始向量x(0)(3,2, 1)T.2 4 4 x123 3 3

13、 x234 42x32解：>>edit Jacobi.mfunctionx iter=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=2 4 -4;3 3 3;4 4 2;> > b=2 -3 -2'> >x0=3 2 -1'> > x,iter=Jacobi(A,x0,b,1e-4)

14、 x =1-1-1 iter =38.用牛顿法求解方程 xcosx 2 0在Xo 2附近的根.解：> > edit Newton.mfunction x iter=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; end end> > edit ksf6.mfunction v=ksf6(x)v=x*cos(x)+2;> > edit ksg6.mfunction z=k

15、sg(y) z=y.A5+y-1;> > x iter=Newton('ksf6','ksg6',2,1e-4) x =2.4988iter =3按模最小9 .分别用改进乘哥法、反哥法计算矩阵A的按模最大特征值及其对应的特征向量 特征值及其对应的特征向量.1266A61626216解：> > edit ep.mfunction t,x=ep(A,x0,tol) tv0 ti0=max(abs(x0); lam0=x0(ti0);x0=x0./lam0;x1=A*x0;tv1 ti1=max(abs(x1);lam1=x1(ti1);x1=

16、x1./lam1;while(abs(lam0-lam1)>tol) x0=x1;lam0=lam1;x1=A*x0;tv1 ti1=max(abs(x1);lam1=x1(ti1);x1=x1./lam1;endt=lam1;x=x1;> > A=12 6 -6;6 16 2;-6 2 16;> >x0=1 0.5 -0.5'> >tol=1e-4;>>t,x=ep(A,x0,tol)t =21.5440x =1.00000.7953-0.7953> > edit fanep.m functiont,x=fanep(A

17、,x0,tol) tv0 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)>tol) x0=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 0.5 -0.5'> >tol=

18、1e-4;t,x=fanep(A,x0,tol) t =4.4560x =1.0000-0.62870.628710 .将积分区间n等分,用复合梯形求积公式计算定积分3Jlx2dx,比较不同n值时的误i差(画出平面上10g(n)-10g(Error)图).解：>> edit ksf7.mfunction v=ksf7(x)v=sqrt(1+x.A2);>> I=quad('ksf7',1,3);>> n=1:100;>>fori=1:100x(i)=tixingjifen('ksf7',1,3,i);error(i

19、)=abs(I-x(i);end>>p1ot(1og10(n),1og10(error(n)11.用SOR迭代方法对下面方程组求解,取初始向量x(0)(1,0, 1)TO21 0 x11131 x2801 2x35解：> > edit sor.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-om

20、ega)*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-4) x = 1.9999 3.0000 -1.0000 iter = 5function splx=-1 0 1 2y=0 -1 0 1 0 -1pp=csape(x,y：complete)breaks,coefs,npolys,ncoefs,dim=unmkpp(pp) xh=-1:0.1:2if -1<=xh

21、<=0yh=coefs(1,1)*(xh+1).A3+coefs(1,2)*(xh+1).A2+coefs(1,3)*(xh+1)+coefs(1,4)elseif 0<=xh<=1yh=coefs(2,1)*(xh).A3+coefs(2,2)*(xh).A2+coefs(2,3)*(xh)+coefs(2,4)elseif 1<=xh<=2yh=coefs(3,1)*(xh-1)A3+coefs(3,2)*(xh-1)A2+coefs(3,3)*(xh-1)+coefs(3,4) elsereturnendplot(xh,yh,'r+')13.

22、二分法程序function x=bisect(f,a,b,tol)if nargin<4tol=1e-12endfa=feval(f,a)fb=feval(f,b)while abs(a-b)>tolx=(a+b)/2 fx=feval(f,x)if sign(fx)=sign(fa)a=xfa=fxelseif sign(fx)=sign(fb) b=xfb=fxelsereturnendend其他>> A=1 2 3 4;5 6 7 8;>> AA =12345678> > m n=size(A)m =2n =4> > x=1 2

23、 3 4 5;>>x x =12345>>length(x) ans =5> > A=2 3 4;1 1 9; 1 2 -6;> > AA =23411912-6>> L U=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>