数值分析重要算法的matlab程序

数值分析重要算法的matlab程序

插值多项式:拉格朗日插值

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

插值多项式:牛顿插值（以数值实验3.2）为例

function y=ex32 n = 21;

x = linspace(-5,5,n)'; h = (5-(-5))/(n-1); y = 1./(1+x.^2);

% form the differences table for j = 2:n,

y(1:n+1-j,j) = diff(y(1:n+2-j,j-1))./(x(j:n)-x(1:n+1-j)); end

% newton coeff y = y(1,:); pz = [ ];

v = linspace(-5,5,80); for t = v,

z = y(n);

for j = n-1:-1:1,

z = z * (t-x(j))+y(j); end

pz=[pz z]; end

plot(v,pz,'r+-',v,1./(1+v.^2),'g--');

数值积分：梯形求积公式求积分

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)+sum(y(2:n))+0.5*y(n+1)); 数值积分：抛物型求积公式求积分 function I=fsimpsion(fun,a,b,n) h=(b-a)/n;

x=linspace(a,b,2*n+1); y=feval(fun,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));

追赶法解三对角方程组

function x=tridiagsolver(a,b) [n,n]=size(a); l=zeros(1,n); y=zeros(1,n); u=zeros(1,n); 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 (i

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 [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:1000;

x=(D-L)\\(b+U*x);

error=norm(b-A*x)/norm(b); if(error

牛顿法求方程的根

function [x, it, convg]=newton(x0, f, g, maxit, tol) %find the zero of function f, with gradient g provided %Usage: [x, it, convg] = newton(x0, f, g, maxit, tol) if nargin<5, tol = 1e-5; if nargin<4, maxit=100; end end x=x0;

fx=feval(f,x); convg = 0; it=1;

while ~convg, it=it+1;

if norm(fx)

fprintf('Newton Iteration successes!!\\n'); convg=1; return; end

d=-feval(g,x)\\fx; lambda=1; lsdone=0; while ~lsdone,

xn=x+lambda*d'; fn=feval(f,xn); if abs(fn)

lambda=1/2*lambda; if lambda<=eps, convg=-1;

error('line search fails!!'); end end end x=xn; fx=fn;

if it > maxit, convg=0;

error('Newton method needs more iterations!!'); end end

龙格-库塔方法求解微分方程数值解

function [x,y]=rk4(ef,tspan,y0,n) y=zeros(1,n+1); y(1)=y0; a=tspan(1); b=tspan(2); h=(b-a)/n; x=a:h:b;

c1=[1 2 2 1]'/6; for i=1:n,

k(1)=h*feval(ef,x(i),y(i));

k(2)=h*feval(ef,x(i)+0.5*h,y(i)+0.5*k(1)); k(3)=h*feval(ef,x(i)+0.5*h,y(i)+0.5*k(2)); k(4)=h*feval(ef,x(i)+h,y(i)+k(3)); y(i+1)=y(:,i)+k*c1; end

乘幂法求特征值和特征向量

function [t,y]=eigIPower(a,xinit,ep) v0=xinit;

[tv,ti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0;

flag=0;

while (flag==0) v1=a*u0;

[tv,ti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1;

err=abs(lam0-lam1); if(err<=ep) flag=1; end

lam0=lam1; end t=lam1; y=u0;

反幂法求特征值和特征向量

function [t,y]=inverPower(a,xinit,ep) v0=xinit;

[tv,ti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0;

while (flag==0) [L,U]=lu(a); yy=L\%u0; v1=U\\yy;

[tv,ti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1;

err=abs(lam0^-1-lam1^-1); if(err<=ep) flag=1; end

lam0=lam1; end

t=lam1^-1; y=u0;

结合原点平移的反幂法求特征值和特征向量

function [t,y]=inverPowerM(a,xinit,p,ep) n=size(a); v0=xinit;

[tv,ti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0;

while (flag==0)

[L,U]=lu(a-p*eye(n)); yy=L\%u0; v1=U\\yy;

[tv,ti]=max(abs(v1));

lam1=v1(ti); u0=v1/lam1;

err=abs(lam0^-1-lam1^-1); if(err<=ep) flag=1; end

lam0=lam1; end

t=lam1^-1; y=u0;