数值分析各种代码.doc

追赶法function x=tridiagsolver(a,b) n,n=size(a); for i=1:n if i=1 l(i)=a(i,i); y(i)=b(i)/l(i); else l(i)=a(i,i)-a(i,i-1)*u(i-1); y(i)=(b(i)-y(i-1)*a(i,i-1)/l(i); end if in u(i)=a(i,i+1)/l(i); end end x(n)=y(n); for j=n-1:-1:1 x(j)=y(j)-u(j)*x(j+1); end 拉格朗日function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:n xp=x(1:i-1 i+1:n); yh=yh+y(i)*prod(c1*xh-xp*c2)./(x(i)-xp*c2);end线性x=x1,x2 y=y1.y2 xh=xh抛物线X=x1,x2,x3 y=y1，y2，y3 xh=xh牛顿差商（输入x,y为列向量）function p,q=chashang(x,y)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1 p(1:n+2-j,j)=diff(p(1:n+3-j,j-1)./(x(j-1:n)-x(1:n+2-j);endq=p(1,2:n+1);三次样条x=0 1 2 3;y=0.2 0 0.5 2.0 1.5 -1;pp=csape(x,y,complete)breaks,coefs,npolys,ncoefs,dim=unmkpp(pp)最小二乘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.1 -0.45 0.27 0.1 -0.29 0.24 0.56 1;A=log(x) cos(x) exp(x);Z=Ay;a0=Z(1)a1=Z(2)a2=Z(3)x=0 0.25 0.50 0.75 1.00;y=1.00 1.284 1.6487 2.1170 2.7183;p=polyfit(x,y,2)a2=p(1)a1=p(2)a0=p(3)复合中点function I=fmid(fun,a,b,n)h=(b-a)/n;x=linspace(a+h/2,b-h/2,n);y=feval(fun,x);I=h*sum(y);复合梯形function I=ftrapz(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,n+1);y=feval(fun,x);I=h*(0.5*y(1)+0.5*y(n+1)+sum(y(2:n);复合辛普森function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(fun,x);I=h/6*(y(1)+y(2*n+1)+2*sum(y(3:2:2*n-1)+4*sum(y(2:2:2*n);雅克比迭代function x,iter=jacobi(A,b,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for iter=1:500 x=D(b+U*x+L*x); error=norm(b-A*x)/norm(b); if (errortol) break; endendGS迭代function x,iter=GS(A,b,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for iter=1:500 x=(D-L)(b+U*x); error=norm(b-A*x)/norm(b); if (errortol) break; endendSOR迭代function x,iter=SOR(A,b,omega,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for iter=1:500 x=(D-omega*L)(omega*b+(1-omega)*D*x+omega*U*x); error=norm(b-A*x)/norm(b); if (errortol) break; endend二分法function v=f(x)v=x3-x-1;function x=erfenfa(f,a,b,tol)if nargintol x=(a+b)/2; fx=feval(f,x); if sign(fx)=sign(fa) a=x; fa=fx; else if sign(fx)=sign(fb) b=x; fb=fx; end endend x=erfenfa(f,1,2)牛顿法function x,it=newton(f,g,x0,tol)it=0;done=0;while done x=x0-feval(f,x0)/feval(g,x0); it=it+1; done=(norm(x-x0)=tol); if done x0=x; endend function r=f(x)r=polyval(1,2,10,-20,x); function r=g(x)p=polyder(1,2,10,-20);r=polyval(p,x);牛顿下山法乘幂法function t,y=chengmifa(a,xinit,ep)v0=xinit;tv,ti=max(abs(v0);lam0=v0(ti);u0=v0/lam0;flag=0;while flag v1=a*u0; tv,ti=max(abs(v1); lam1=v1(ti); u0=v1/lam1; err=abs(lam0-lam1); if errep flag=1 end lam0=lam1;endt=lam1;y=u0; 反幂法function t,y=fanmifa(a,xinit,ep)v0=xinit;tv,ti=max(abs(v0);lam0=v0(ti);u0=v0/lam0;flag=0;while flag v1=inv(a)*u0; tv,ti=max(abs(v1); lam1=v1(ti); u0=v1/lam1; err=abs(lam0-lam1); if errep flag=1; end lam0=lam1;endt=1/lam1;y=u0;改进欧拉法function x,y=odeIEuler(f,y0,a,b,