各种迭代法编程

雅可比迭代法：

function x=jacobi(a,b,p,delta,n)

%a为n维非奇异矩阵；b为n维值向量

%p为初值；delta为误差界；n为给定的迭代最高次数

N=length（b）;

for k=1:n

for j=1：N

x(j)=(b(j)-a(j,[1:j-1,j+1:N])*p([1:j-1,j+1:N]))/a(j,j);

end

err=abs(norm(x’-p));

p=x’;

if(err

break;

end

end

p %显示迭代过程

x=x’；

k,err

高斯塞德尔法迭代：

function x=saidel(a,b,p,delta,n)

%a为n维非奇异矩阵；b为n维值向量

%p为初值；delta为误差界；n为给定的迭代最高次数

N=length（b）;

for k=1:n

for j=1：N

if j==1

x(1)=(b(1)-a(1,2:N)*p(2:N))/a(1,1);

else if j=N

x(N)=(b(N)-a(N,1:N-1)*(x(1:N-1))’)/a(N,N);

else

x(j)=(b(j)-a(j,(1:j-1)*x(1:j-1）-a(j,j+1:N)*p(j+1:N))/a(j,j);

end

end

err=abs(norm(x’-p));

p=x’;

if(err

break;

end

end

x=x’；

k,err

不动点迭代法：

function [x,k,err,p]=ddf(f,x0,tol,n)

%ddl.m为用迭代法求非线性方程的解

%f为给定的迭代函数；x0为给定的初始值

%tol为给定的误差界；n为所允许的最大迭代次数

%k为迭代次数；x为不动点的近似值；err为误差

p(1)=x0;

for k=2:n

p(k)=feval(f,p(k-1));

k,

err=abs(p(k)-p(k-1))

x=p(k);

if(err

break;

end

if k==n

disp('迭代超过最大次数！')

end

end

x=p'

牛顿法：

function [x,k,err,y]=Newtun(f,df,x0,tol,n)

%Newtun.m为用迭代法求非线性方程的解

%f为给定的非线性方程；df为f的微分方程；x0为给定的初始值%tol为给定的误差界；n为所允许的最大迭代次数

%k为迭代次数；x为不动点的近似值；err为误差

%x为牛顿迭代法得到得近似解

y(1)=x0;

for k=1:n

x=x0-feval(f,x0)/feval(df,x0);

err=abs(x-x0);

x0=x;

if(err

break;

end

end

必要编辑M文件qfun.m,代码如下：

function y=qfun(x);

y=x^3-3*x-1;

弦截法：

function [x,err,k,y]=xjf(f,x0,x1,tol,n)

%xjf.m为用弦截法迭代法求非线性方程的解

%f为给定的非线性方程；x0,x1为给定的初始值

%tol为给定的误差界；n为所允许的最大迭代次数

%k为迭代次数；x为牛顿迭代法的近似值；err为x1-x0的绝对值

y(1)=x0;

for k=1:n

x=x1-feval('li6_5fun',x1)*(x1-x0)/(feval('li6_5fun',x1)-feval('li6_5fun',x0));

err=abs(x-x1);

x0=x1;

x1=x;

if(err

break;

end

end

必要编辑M文件li6_5.m,代码如下：

function y=li6_5(x);

y=x^3-3*x-1;

复化梯形公式matlab:

function t=tixing(f,a,b,n)

h=(b-a)/n;

sum=0;

for k=0:n-1

x=a+k*h;

sum=sum+feval(f,x);

end

t=h/2*(feval(f,a)+feval(f,b)+2*sum);

运行程序结果：

>> format long

>> tixing(inline('x./(x.^2+4)'),0,1,8)

ans =

0.111402354529548

复化辛普森公式matlab:

function y=xinpusen(f,a,b,n)

h=(b-a)/n;

sn=h/6*(feval(f,a)+feval(f,b));

s=0;

for i=0:n-1

s=s+4*feval(f,a+h/2+i*h)+2*feval(f,a+i*h);

end

y=sn+h*s/6;

运行结果：

>> xinpusen(inline('x./(x.^2+4)'),0,1,8)

ans =