各种迭代法编程

雅可比迭代法：

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

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)

Yl.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

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

必要编辑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

必要编辑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 =

0.111571813252631

改进欧拉法:

function y=gjEuler(f,x0,y0,xn,N)

%gjEuler.m函数为改进的欧拉法求微分方程的解 %f为一阶常微分方程的一般表达式的右端函数 %x0,y0为初始条件

%xn为取值范围的一个端点 %N为区间的个数

%y为求解微分方程组的解 x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;

h=(xn-x0)/N; for n=1:N

x(n+1)=x(n)+h;

z0=y(n)+h*feval(f,x(n),y(n));

y(n+1)=y(n)+h/2*(feval(f,x(n),y(n))+feval(f,x(n+1),z0)); end

T=[x',y']

梯形法：

function y=TXF(f,x0,y0,xn,N)

%TXF.m函数为改进的梯形法求微分方程的解 %f为一阶常微分方程的一般表达式的右端函数 %x0,y0为初始条件

%xn为取值范围的一个端点 %N为区间的个数

%y为求解微分方程组的解 x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;

h=(xn-x0)/N; for n=1:N

x(n+1)=x(n)+h;

y(n+1)=y(n)+h/2*(feval(f,x(n),y(n))+feval(f,x(n+1),y(n+1))); end T=[x',y']

function y=LG4(f,x0,y0,xn,N) %LG4.m函数为四阶龙格库塔法求微分方程的解 %f为一阶微分方程

%x0,y0为左右端点 %xn为给定的初始值 %N为给定的迭代步长 %y为求解微分方程组的解 x=zeros(1,N+1); y=zeros(1,N+1); h=(y0-x0)/N; T=x0:h:y0; y(1)=xn;

for n=1:N

k1=h*feval(f,T(n),y(n));

k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+h/2,y(n)+k2/2); k4=h*feval(f,T(n)+h,y(n)+k3);

y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6; end

R=[T',y']

直接三角分解：

function dxl=dxllu(a,b) format rat

[n n] =size(a); RA=rank(a); if RA~=n

disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:') RA

hl=det(a) return end

if RA==n for p=1:n

h(p)=det(a(1:p, 1:p)); end

hl=h(1:n); for i=1:n

if h(i)==0

disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：') hl RA return end end

本文来源：https://www.bwwdw.com/article/c202.html