数值分析重要算法的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<n)
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<tol)
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)<tol,
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)<abs(fx),
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=[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;。