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

1、同济大学数值分析matlab编程题汇编同济大学数值分析matlab编程题汇编MATLAB 编程题库1.下面的数据表近似地满足函数，请适当变换成为线性最小二乘问题，编程求最好的系数，并在同一个图上画出所有数据和函数图像.解：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;yh=(a.*xh+b)./(1+c.*xh.

2、2);plot(x,y,r+,xh,yh,b*)2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数精度为的近似根，并写出调用方式：文件一文件二function v = f(x) v = x .* log(x) - 1;function z = g(y) z = y.5 + y - 1;解： edit gexianfa.mfunction 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)-fe

3、val(f,x0);x0=x1;x1=x;end edit f.mfunction v=f(x)v=x.*log(x)-1; edit g.mfunction z=g(y)z=y.5+y-1; x1 iter1=gexianfa(f,1,3,1e-10)x1 = 1.7632iter1 = 6 x2 iter2=gexianfa(g,0,1,1e-10)x2 = 0.7549iter2 = 83.使用GS迭代求解下述线性代数方程组：解： edit gsdiedai.mfunction x iter=gsdiedai(A,x0,b,tol)D=diag(diag(A);L=D-tril(A);U

4、=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=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.0910 1.2372 0.9802iteriter = 64.用四阶Range-kutta方法求解下述常微分方程初值问题（取步长h=0.01）解： edit ksf2.mfunction v=ksf2(x,y)v=y+exp(x)+x.

5、*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+k4)./6;endy调用函数方法 edit Rangekutta.mfunction x y=Rangekutta(f,a,b,h,y0)x=a:h:b;n=(b

6、-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)./6;end x y=Rangekutta(ksf2,1,2,0.01,2);y5.取，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在上求解初值问题。解： edit Euler.m

7、function 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.mfunctionx 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)./2;end edit ksf3.mfu

8、nction v=ksf3(x,y)v=x.3-y./x;x y=Euler(ksf3,1,2,0.2,0.4)x = 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.4)x = 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程序，计算下面积分的近

9、似值，区间等分。编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分。、解： edit tixingjifen.mfunction 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*

10、sum(y(3:2:2*n-1)+4*sum(y(2:2:2*n)+y(2*n+1); edit ksf4.mfunction v=ksf4(x)v=1./(x.2+1);tixingjifen(ksf4,0,1,20)ans = 0.7853simpson(ksf4,0,1,10)ans = 0.7854 edit ksf5.mfunction v=ksf5(x)if(x=0) v=1;else v=sin(x)./x;end(第二个函数ksf5调用求积函数时，总显示有错误：“NaN”，还没调试好。见谅！)7.用迭代方法对下面方程组求解，取初始向量。解：edit Jacobi.mfuncti

11、onx 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)x = 1 -1 -1iter = 38.用牛顿法求解方程在附近的根。解： edit Newton.mfunction x iter=Newton(f,g,x0,tol)it

12、er=0;done=0while donex=x0-feval(f,x0)/feval(g,x0);done=norm(x-x0) edit ksf6.mfunction v=ksf6(x)v=x*cos(x)+2; edit ksg6.mfunction z=ksg(y)z=y.5+y-1; x iter=Newton(ksf6,ksg6,2,1e-4)x = 2.4988iter = 39.分别用改进乘幂法、反幂法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。解： edit ep.mfunction t,x=ep(A,x0,tol)tv0 ti0=max(

13、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)tol)x0=x1; lam0=lam1;x1=A*x0; tv1 ti1=max(abs(x1); lam1=x1(ti1);x1=x1./lam1;endt=lam1;x=x1; edit fanep.mfunctiont,x=fanep(A,x0,tol)tv0 ti0=max(abs(x0);lam0=x0(ti0);x0=x0./lam0;x1=Ax0;tv1 ti1=max

14、(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=1e-4;t,x=ep(A,x0,tol)t = 21.5440x = 1.0000 0.7953 -0.7953 A=12 6 -6;6 16 2;-6 2 16;x0=1 0.5 -0.5;tol=1e-4;

15、t,x=fanep(A,x0,tol)t = 4.4560x = 1.0000 -0.6287 0.628710.将积分区间n等分，用复合梯形求积公式计算定积分,比较不同值时的误差（画出平面上 log(n)-log(Error)图）解： edit ksf7.mfunction v=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);endplot(log10(n),log10(error(n)11.用迭代方法对下面方程组求解，取初始向

16、量。解： 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-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-4)x = 1.9999 3.0000 -1.0000it

17、er = 512.编程实现求解满足下列条件的区间-1,2上的三次样条函数S(x),并画出此样条函数的图形:xi -1 0 1 2f(xi) -1 0 1 0f(xi) 0 -1function 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=0 yh=coefs(1,1)*(xh+1).3+coefs(1,2)*(xh+1).2+coefs(1,3)*(xh+1)+coefs(1,4)elseif 0=xh=1 yh=co

18、efs(2,1)*(xh).3+coefs(2,2)*(xh).2+coefs(2,3)*(xh)+coefs(2,4)elseif 1=xh=2 yh=coefs(3,1)*(xh-1).3+coefs(3,2)*(xh-1).2+coefs(3,3)*(xh-1)+coefs(3,4)else returnendplot(xh,yh,r+)13.二分法程序function x=bisect(f,a,b,tol)if nargintol x=(a+b)/2 fx=feval(f,x) if sign(fx)=sign(fa) a=x fa=fx elseif sign(fx)=sign(fb

19、) b=x fb=fx else return endend其他 A=1 2 3 4;5 6 7 8; AA = 1 2 3 4 5 6 7 8 m n=size(A)m = 2n = 4 x=1 2 3 4 5;xx = 1 2 3 4 5length(x)ans = 5 A=2 3 4;1 1 9; 1 2 -6; AA = 2 3 4 1 1 9 1 2 -6 L U=lu(A); LL = 1.0000 0 0 0.5000 1.0000 0 0.5000 -1.0000 1.0000 UU = 2.0000 3.0000 4.0000 0 -0.5000 7.0000 0 0 -1.

20、0000 x=linspace(0,1,11);xans = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.90001.0000 x=1 2 3 4; y=6 11 18 27; p=polyfit(x,y,2)p =1.0000 2.0000 3.0000diag(ones(4,1),1)ans = 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 012.编程实现求解满足下列条件的区间-1,2上的三次样条函数S(x),并画出此样条函数的图形:xi -1 0 1 2f(xi

21、) -1 0 1 0f(xi) 0 -1function 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=0 yh=coefs(1,1)*(xh+1).3+coefs(1,2)*(xh+1).2+coefs(1,3)*(xh+1)+coefs(1,4)elseif 0=xh=1 yh=coefs(2,1)*(xh).3+coefs(2,2)*(xh).2+coefs(2,3)*(xh)+coefs(2,4)elseif 1=xh=2 yh=coefs(3,1)*(xh-1).3+coefs(3,2)*(xh-1).2+coefs(3,3)*(xh-1)+coefs(3,4)else returnendplot(xh,yh,r+)13.二分法程序function x=bisect(f,a,b,tol)if nargintol x=(a+b)/2 fx=feval(f,x) if sign(fx)=sign(fa) a=x fa=fx elseif sign(fx)=sign(fb) b=x fb=fx else return endend