Matlab经典程序 现代数值计算

追赶法程序：

function x=zgf(A,b)

[n,n]=size(A);

x=zeros(n,1);

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)-A(i,i-1)*y(i-1))/l(i);

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)-y(j+1)*u(j);

end

拉格朗日插值程序：

function yh=lage(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]);

B=(c1*xh-xp'*c2)./(x(i)-xp'*c2)

yh=yh+y(i)*prod(B);

end

三次样条函数：

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) 复合梯形公式程序：

function I=ti(f,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=simp(f,a,b,n)

h=(b-a)/n;

x=linspace(a,b,2*n+1);

y=feval(f,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,iter]=jacobi(A,b,tol)

D=diag(diag(A))

L=-tril(A,-1)

U=-triu(A,1)

x=zeros(size(b));

for iter=1:500

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

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

if(error

break;

end

end

GS迭代程序：

function [x,iter]=gs(A,b,tol)

D=diag(diag(A))

L=-tril(A,-1)

U=-triu(A,1)

x=zeros(size(b));

for iter=1:500

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

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

if(error

break;

end

end

SOR迭代程序：

function [x,iter]=sor(A,b,w,tol)

D=diag(diag(A))

L=-tril(A,-1)

U=-triu(A,1)

x=zeros(size(b));

for iter=1:500

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

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

if(error

break;

end

end

牛顿法求解非线性方程的根：function f=f(x)

f=x.^2+sin(10*x)-1;

function df=df(x)

g=2*x+10*cos(10*x);

牛顿法程序开始：

function x=newton(f,g,a,b,tol)

x0=(a+b)/2;

flag=0;

whie(flag==0)

x1=x0-feval(f,x0)/feval(g,x0);

err=abs(x1-x0);

if(err

flag=1

end

x0=x1

end

x=x0; 割线法：

function f=f(x)

f=x.^2+sin(10*x)-1;

function x=gexian(f,a,b,tol,)

x0=(a+b)/2;

x1=(a+b)/2+1;

flag=0;

while(flag==0)

x2=x1-feval(f,x1)*(x1-x0)/(feval(f,x1)-

feval(f,x0));

err=abs(x2-x1);

if(err

flag=1

end

x1=x2

end

x=x1

改进乘幂法程序：

function [t,y]=chengmi(a,x0,ep)

v0=x0;

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