同济大学数值分析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=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工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数精度为 的近似根，并写出调用方式： 文件一 文件二 function v = f(x) v = x .* log(x) - 1; function z = g(y) z = y.^5 + y - 1;     解： >> edit gexianfa.m function [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.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8 3.使用GS迭代求解下述线性代数方程组： 解： >> edit gsdiedai.m function [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=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); >>x x = -3.0910 1.2372 0.9802 >>iter iter = 6 4.用四阶Range-kutta求解下述常微分方程初值问题（取步长h=0.01） 解： >> edit 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+k4)./6; end >>y 调用函数方法 >> edit Rangekutta.m function [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)./6; end >> [x y]=Rangekutta('ksf2',1,2,0.01,2); >>y 5.取 ，请编写Matlab程序，分别用欧拉方法、改进欧拉方法在 上求解初值问题。 解： >> edit Euler.m 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.m function[x 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.m function 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.0000 y = 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.0000 y = 0.4000    0.5895    0.9278    1.4615    2.2464    3.3466 6.请编写复合梯形积分公式的Matlab程序，计算下面积分的近似值，区间等分 。 编写辛普森积分公式的Matlab程序，计算下面积分的近似值，区间等分 。 、 解： >> edit tixingjifen.m function 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.m function 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.m function v=ksf4(x) v=1./(x.^2+1); >>tixingjifen('ksf4',0,1,20) ans = 0.7853 >>simpson('ksf4',0,1,10) ans = 0.7854 >> edit ksf5.m function v=ksf5(x) if(x==0) v=1; else v=sin(x)./x; end (第二个函数‘ksf5’调用求积函数时，总显示有错误：“NaN”，还没调试好。见谅！) 7.用 迭代方法对下面方程组求解，取初始向量 。 解： >>edit Jacobi.m function[x 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 -1 iter = 3 8.用牛顿法求解方程 在 附近的根。 解： >> edit Newton.m function [x iter]=Newton(f,g,x0,tol) iter=0; done=0 while ~done x=x0-feval(f,x0)/feval(g,x0); done=norm(x-x0)<=tol; iter=iter+1; if  ~done,x0=x; end end >> edit ksf6.m function v=ksf6(x) v=x*cos(x)+2; >> edit ksg6.m function z=ksg(y) z=y.^5+y-1; >> [x iter]=Newton('ksf6','ksg6',2,1e-4) x = 2.4988 iter = 3 9.分别用改进乘幂法、反幂法计算矩阵A的按模最大特征值及其对应的特征向量、按模最小特征值及其对应的特征向量。 解： >> edit ep.m function [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));