数值分析重要算法的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-1i+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

break;

end

end

牛顿法求方程的根

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)

lsdone=1;

else 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=[1221]'/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;