同济大学数值分析历年程序解答.docx

第二章线性方程组直接解法追赶法function x=zhuigan(A,d)n n=size(A);for i=1:n if i=1 %这里要注意= l(i)=A(i,i); y(i)=d(i)/l(i); u(i)=A(i,i+1)/l(i); else l(i)=A(i,i)-A(i,i-1)*u(i-1); y(i)=(d(i)-y(i-1)*A(i,i-1)/l(i);ifinu(i)=A(i,i+1)/l(i);endendendx(n)=y(n);fori=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);end %ai,bi,ci在程序里需要按照其在A矩阵中的位置写出，所以要把矩阵写出来，看其位置写相应的项第三章函数插值y1=interp1(x,y,x1,linear) %线性插值，根据x,y表示的函数，计算x1处的函数值y1=interp1(x,y,x1,cubic)%三次多项式插值y1=interp1(x,y,x1,spline)%三次样条插值第四章函数逼近/拟合x=0.24 0.65 0.95 1.24 1.73 2.01 2.23 2.52 2.77 2.99; %任意选取拟合函数的最小二乘程序y=0.23 -0.26 -1.10 -0.45 0.27 0.1 -0.29 0.24 0.56 1;A=log(x) cos(x) exp(x); z=Ay%A矩阵为选取函数各点值构成的矩阵，z为系数矩阵，求解系数，等价于解方程A*z=yx=0 0.25 0.5 0.75 1; %多项式拟合的最小二乘程序y=1 1.284 1.6487 2.117 2.7183;x1=ones(5,1);A=x1 x x.2; z=Ayp=polyfit(x,y,n); %多项式最小二乘拟合dp=polyder(p); %多项式求导y=polyval(dp,x) %多项式求值函数，这里求的是点x处导数值第五章求积公式function I=mid(f,a,b,n) %中点公式 h=(b-a)/n; x=a+h/2:h:b-h/2; %中点自变量的值 y=feval(f,x); I=h*sum(y);function I=trapeze(f,a,b,n) %梯型公式 h=(b-a)/n; x=a:h:b; y=feval(f,x); I=h/2*(y(1)+2*sum(y(2:n)+y(n+1); %f=h/2*f(a)+2f(xi)+f(b)function I=simpson(f,a,b,n) %simpson公式 h=(b-a)/(2*n); x=a:h:b; y=feval(f,x); I=h/3*(y(1)+4*sum(y(2:2:2*n)+2*sum(y(3:2:2*n-1)+y(2*n+1); %f=h/3*f(a)+4f(ou)+2f(ji)+f(b)第六章线性方程组的迭代解Jacobi迭代function x,iter=jacobi(A,b,tol); % Ax=b ,iter迭代次数D=diag(diag(A); % (D-L-U)x=bL=D-tril(A); % Dx=Lx+Ux+bU=D-triu(A); x=zeros(size(b);foriter=1:50000; x=D(L+U)*x+Db; error=norm(A*x-b)/norm(b);if errortolbreak;endend高斯-赛德尔迭代function x,iter=gs(A,b,tol); % Ax=b ,iter迭代次数D=diag(diag(A); % (D-L-U)x=bL=D-tril(A); % (D-L)x=+Ux+bU=D-triu(A);x=zeros(size(b);foriter=1:50000; x=(D-L)U*x+(D-L)b;error=norm(A*x-b)/norm(b);if errortolbreak;endend第七章非线性方程求根二分法function x=erfen(f,a,b,tol) %二分法求解方程f(x)=0，a,b,为给定的初始区间，应满足f(a)*f(b)0 while 1 c=(a+b)/2; if norm(a-b)=tol x=a; return end y1=feval(f,a); y2=feval(f,c); if sign(y1)=sign(y2) a=c; else b=c; end end牛顿法function x,it=Newton(f,df,x0,maxit,tol) %求解方程f(x)=0% it迭代次数, maxit,最大迭代次数, tol,迭代精度 f求解的函数df,f的导数it=0;while 1; x=x0-feval(f,x0)/feval(df,x0); %注意除号方向it=it+1; if norm(x-x0)=maxiterror(Newton method needs more itrations)returnendend第八章矩阵特征值和特征向量幂法function lam x=eigIPower(A,x0,tol) v0=x0; m n=max(abs(x0); lam0=v0(n); u0=v0/lam0;while 1 v1=A*u0; m n=max(abs(v1); lam1=v1(n); u1=v1/lam1;if abs(lam1-lam0)=tollam=lam1; x=u1;returnend u0=u1; lam0=lam1;end%幂法求解Ax=lam*x最大特征值和特征向量 lam最大特征值，x对应的特征向量,x0迭代初始向量，tol迭代精度第九章常微分方程初值问题Euler法function x y=odeEuler(f,y0,a,b,n) y(1)=y0; h=(b-a)/n; x=a:h:b;fori=1:n y(i+1)=y(i)+h*feval(f,x(i),y(i);end% x,自变量值 y,函数值 f导数/增量函数 y0初值a,b,n求解区间开始点，终止点，划分区间数改进Euler法function x y=odeIEuler(f,y0,a,b,n) y(1)=y0; h=(b-a)/n; x=a:h:b;fori=1:n y1=y(i)+h*feval(f,x(i),y(i); y(i+1)=y(i)+0.5*h*(feval(f,x(i),y(i)+feval(f,x(i+1),y1); end%x,自变量值 y,函数值 f导数/增量函数 y0初值a,b,n求解区间开始点，终止点，划分区间数%两点导取平均，用下一个点的预估值算下一个点的导数值Runge-Kutta法function x y=RungeKutta(f,y0,a,b,n)y(1)=y0; h=(b-a)/n; x=a:h:b;fori=1:n y1=h*feval(f,x(i),y(i); y2=h*feval(f,x(i)+0.5*h,y(i)+0.5*y1); y3=h*fev