【免费下载】数值分析的几个简单算法实现

数值分析的几个简单算法实现 Matlab6.5

function x=Aitken(fname,x0,N,e)

k=1;e=1e-5;x2=x0+2*e;

while abs(x2-x0)>e&k

x1=feval(x0);x2=feval(x1);

x2=x2-(x2-x1)^2/(x2-2*x1+x0);

k=k+1;

x2=x0;end;

disp(x2);

if k==N,warning('迭代失败');end;

end

%迭代法diedai.m文件

function x=diedai(fname,x0,e,N)

if nargin<5;N<500;end;

if nargin<4;e=1e-4;end;

x=x0;x0=x+2*e;e=1e-5;k=1;

while abs(x0-x)>e&k

x=feval(fname,x0);

k=k+1;

x0=x;

disp(x);

end

if k==N,warning('迭代失败');end

%二分法erfenfa.m

function x=erfenfa(fname,a,b,e)

y0=feval(fname,a);e=1e-5;

while (b-a)>e

x=(a+b)/2;

y=feval(fname,x);

if y*y0>0

a=x;

else

b=x;

end

disp(x);

end

%牛顿法

function x=newton(fname,dfname,x0,e,N)

if nargin<5;N<500;end;

if nargin<4;e=1e-4;end;

x=x0;x0=x+2*e;k=0;

while abs(x0-x)>e&k

k=k+1;

x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);

disp(x);

end

if k==N,warning('已达迭代次数上限');end

% format long; newton(inline('exp(x)-4*cos(x)'),inline('exp(x)+4*sin(x)'),pi/4,1e-4,20) function t=Rombeg(fname,a,b,e)

%用途：龙贝格法求函数积分

%格式：t=Rombeg(fname,a,b,e) fname是被积函数，a,b分别为下上限，e为精度（默认1e-4）

if nargin<4,e=1e-4;end;

i=1;j=1;h=b-a;

T(i,1)=h/2*(feval(fname,a)+feval(fname,b));

T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:1-h/2+0.001*h))*h/2;

T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);

while abs(T(i+1,i+1)-T(i,i))>e

i=i+1;h=h/2;

T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:1-h/2+0.001*h))*h/2;

for j=1:i

T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);

end

end

T

t=T(i+1,j+1);

% format long;Rombeg(inline('sin(x)./x'),eps,1,1e-6)

function x=gaussseidel(A,b,x0,e,N)

%用途：用向量形式的gauss-seidel迭代解线性方程组Ax=b

if nargin<5,N=500;end

x=x0;x0=x+2*e;

k=0;A1=tril(A);A2=inv(A1);

while norm(x0-x,inf)>e&k

k=k+1;

x0=x;x=-A2*(A-A1)*x0+A2*b;

disp(x')

end

if k==N,warning('已达到迭代次数上限');

end

>> x=gaussseidel(A,b,[1 1 1]',1e-6)

function x=sor(A,b,omega,x0,e,N)

%用途：用分量形式的sor迭代解线性方程组Ax=b n=length(b);

if nargin<6,N=500;end

x=x0;x0=x+2*e;

k=0;L=tril(A,-1);U=triu(A,1);

while norm(x0-x,inf)>e&k

k=k+1

x0=x;

for i=1:n

x1(i)=(b(i)-L(i,1:i-1)*x(1:i-1,1)-U(i,i+1:n)*x0(i+1:n,1))/A(i,i);

x(i)=(1-omega)*x0(i)+omega*x1(i);

end

disp(x')

end

if k==N,warning('已达到迭代次数上限');

end

>> sor(A,b,1.45,[1 1 1]',1e-6)

function x=jacobi(A,b,x0,e)

%用途：用分量形式的sor迭代解线性方程组Ax=b

n=length(b);

x=x0;x0=x+2*e;

k=0;L=tril(A,-1);U=triu(A,1);

while norm(x0-x,inf)>e,

k=k+1;

x0=x;

for i=1:n

x(i)=(b(i)-L(i,1:i-1)*x0(1:i-1,1)-U(i,i+1:n)*x0(i+1:n,1))/A(i,i);

end

disp(k)