数值方法.docx
《数值方法.docx》由会员分享，可在线阅读，更多相关《数值方法.docx（42页珍藏版）》请在冰点文库上搜索。
数值方法
非线性方程数值方法
1、固定点迭代：
求解方程
的近似值，初始值为
，迭代式为
function[k,x1,err,x]=fixpt（g,x0,tol,maxl）
%Input-gistheiterationfunctioninputasastring'g'
%-x0istheinitialguessforthefixedpoint
%-tolisthetolerance
%-maxlisthemaximumnumberofiterations
%Output-kisthenumberofiterationthatwerecarriedout
%-x1istheapproximationtothefixedpoint
%-erristheerrorintheapproximation
%-xcontainsthesequence{xn}
x
（1）=x0;
fork=2:
maxl
x（k）=feval（g,x（k-1））;
err=abs（x（k）-x（k-1））;
relerr=err/（abs（x（k））+eps）;
x1=x（k）;
if（errend
ifk==maxl
disp（'maximumnumberofiterationsexceeded'）
end
x=x';
例如 g=inline（'exp（-x）'）;
[k,p,err,P]=fixpt（g,0.5,0.01,20）
2、二分法：
求解方程
在区间
内的一个根。
前提条件是
是连续的，且
与
的符号相反。
function[c,err,yc]=bisect（f,a,b,delta）
%Input-fisthefunctioninputasastring'f'
%-aandbaretheleftandrightendpoint
%-deltaisthetolerance
%Output-cisthezero
%-yc=f（c）
%-erristheerrorestimateforc
ya=feval（f,a）,yb=feval（f,b）
ifya*yb>0,break,end
maxl=1+round（（log（b-a）-log（delta））/log
（2））
fork=1:
maxl
c=（a+b）/2;yc=feval（f,c）;
ifyc==0
a=c;b=c;
elseifyb*yc>0
b=c;yb=yc;
else
a=c;ya=yc;
end
ifb-aend
c=（a+b）/2;err=abs（b-a）;yc=feval（f,c）;
例：
f=inline（'x*sin（x）-1'）
[c,err,yc]=bisect（f,0,2,0.01）
3、牛顿－拉夫申迭代：
使用初始近似值
，利用迭代式
，计算函数
的根的近似值。
function[x0,err,k,y]=newton（f,df,x0,delta,epsilon,maxl）
%Input-fistheobjectfunctioninputasastring'f'
%-dfisthederivativeoffinputasastring'df'
%-x0istheinitialapproximationtoazerooff
%-deltaisthetoleranceforx0
%-epsilonisthetoleranceforthefunctionvaluesy
%-maxlisthemaximumnumberofiterations
%Output-xistheNewton-Raphsonapproximationtothezero
%-erristheerrorestimateforx0
%-kisthenumberofiterations
%-yisthefunctionvaluef（x0）
fork=1:
maxl
x1=x0-feval（f,x0）/feval（df,x0）;
err=abs（x1-x0）;
relerr=2*err/（abs（x1）+delta）;
x0=x1;
y=feval（f,x0）;
if（errend
例：
f=inline（'x^2-5'）,df=inline（'2*x'）
[x,err,k,y]=newton（f,df,2,0.001,0.001,20）
4、割线法：
使用初始近似值
和
，利用迭代式
，计算函数
的根的近似值。
function[x1,err,k,y]=secant（f,x0,x1,delta,epsilon,maxl）
%Input-fistheobjectfunctioninputasastring'f'
%-x0andx1istheinitialapproximationstoazero
%-deltaisthetoleranceforx1
%-epsilonisthetoleranceforthefunctionvaluesy
%-maxlisthemaximumnumberofiterations
%Output-x1isthesecantmethodapproximationtothezero
%-erristheerrorestimateforx1
%-kisthenumberofiterations
%-yisthefunctionvaluef（x1）
fork=1:
maxl
x2=x1-feval（f,x1）*（x1-x0）/（feval（f,x1）-feval（f,x0））;
err=abs（x2-x1）;
relerr=2*err/（abs（x1）+delta）;
x0=x1;x1=x2;
y=feval（f,x1）;
if（errend
例：
f=inline（'x^3-3*x+2'）
[x1,err,k,y]=secant（f,-2.6,-2.4,0.01,0.001,20）
5、Steffensen加速法：
给定初始近似值
，快速寻找固定点方程
的解，假设
和
是连续的，而且
，通常的固定点迭代缓慢（线性）收敛到
。
function[P,Q]=steff（f,df,p0,delta,epsilon,maxl）
%Input-fistheobjectfunctioninputasastring'f'
%-dfisthederivativeoffinputasastring'df'
%-p0istheinitialapproximationtoazerooff
%-deltaisthetoleranceforp0
%-epsilonisthetoleranceforthefunctionvaluesy
%-maxlisthemaximumnumberofiterations
%Output-PistheSteffensenapproximationtothezero
%-QisthematrixcontainingtheSteffensensequence
%InitializethematrixR
R=zeros（maxl,3）;
R（1,1）=p0;
fork=1:
maxl
forj=2:
3
%DenominatorinNewton-Raphsonmethodiscalculated
nrodenom=feval（df,R（k,j-1））;
%CalculateNewton-Raphsonapproximations
ifnrodenom==0
'divisionbyzeroinNewton-Raphsonmethod'
break
else
R（k,j）=R（k,j-1）-feval（f,R（k,j-1））/nrodenom;
end
%DenominatorinAitken'sAccelerationprocesscalculated
aadenom=R（k,3）-2*R（k,2）+R（k,1）;
%CalculateAitken'sAccelerationapproximations
ifaadenom==0
'divisionbyzeroinAitken''sAcceleration'
break
else
R（k+1,1）=R（k,1）-（R（k,2）-R（k,1））^2/aadenom;
end
end
%Endprogramifdivisionbyzerooccurred
if（nrodenom==0）|（aadenom==0）
break
end
%Stoppingcriteriaareevluated
err=abs（R（k,1）-R（k+1,1））;
relerr=err/（abs（R（k+1,1））+delta）;
y=feval（f,R（k+1,1））;
if（err%PandthematrixQaredetermined
P=R（k+1,1）;Q=R（1:
k+1,:
）;
break
end
end
end
例：
f=inline（'x^2-5'）,df=inline（'2*x'）
[P,Q]=steff（f,df,2,0.001,0.001,20）
6、Muller法：
给定三个初始值
，求方程
的根。
Muller法是一种迭代方法，构造一条抛物线经过初始三点，然后利用二次函数求下一个近似值的二次根。
可以证明当接近一个单根时，Muller法比割线法收敛得快，基本上与牛顿法一样快。
此方法可用来求函数的实数零点和复数零点，且可用于为复杂的算式编程序。
设
是根的最佳近似值。
改变变量
。
使用差分
和
。
设包含
的二次多项式为：
。
根据三点
可得到三个包含
的三个方程，解得
，并利用定义
和
，可得
（1）
二次式的根为
（2）
为确保方法的稳定性，需选择式是绝对值最小的根。
如果
，使用带正号的根；如果
，使用带负号的根。
则最新的根为
。
为更新迭代，需要从
中选择最靠近
的两点为新的
（即放弃离
最远的一点）。
function[p,y,err]=muller（f,p0,p1,p2,delta,epsilon,maxl）
%Input-fistheobjectfunctioninputasastring'f'
%-p0,p1,andp2aretheinitialapproximations
%-deltaisthetoleranceforp0,p1,andp2
%-epsilonisthetoleranceforthdfunctionvaluesy
%-maxlisthemaximumnumberofiterations
%Output-pistheMullerapproximationtothezerooff
%-yisthefunctionvaluey=f（p）
%-erristheerrorintheapproximationofp
%InitializethematrixesPandY
P=[p0p1p2];
Y=feval（f,P）;
%Calculateaandbinformula
（1）
fork=1:
maxl
h0=P
（1）-P（3）;h1=P
（2）-P（3）;e0=Y
（1）-Y（3）;e1=Y
（2）-Y（3）;c=Y（3）;
denom=h1*h0^2-h0*h1^2;
a=（e0*h1-e1*h0）/denom;
b=（e1*h0^2-e0*h1^2）/denom;
%Suppressanycomplexroots
ifb^2-4*a*c>0
disc=sqrt（b^2-4*a*c）;
else
disc=0;
end
%Findthesmallestroot
（2）
ifb<0
disc=-disc;
end
z=-2*c/（b+disc）;
p=P（3）+z;
%sorttheentriesofPtofindthetwoclosesttop
ifabs（p-P
（2））（1））
Q=[P
（2）P
（1）P（3）];P=Q;Y=feval（f,P）;
end
%ReplacetheentryofPthatwasfarthestfromPwithp
P（3）=p;Y（3）=feval（f,P（3））;
y=Y（3）;
%Determinestoppingcriteria
err=abs（z）;
relerr=err/（abs（p）+delta）;
if（errbreak
end
end
例：
f=inline（'x.^3-3.*x+2'）;[p,y,err]=muller（f,1.4,1.3,1.2,0.01,0.001,20）
线性方程组数值方法
1、回代：
用回代法求解上三角线性方程组
，必须满足系数矩阵的对角元素非零的条件。
functionx=backsub（A,b）
%Input-Aisannxnupper-triangularnonsingularmatrix
%-bisannx1matrix
%Output-xisthesolutiontothelinearsystemAx=b
%Findthedimensionofbandinitializex
n=length（b）;x=zeros（n,1）;
x（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（b（k）-A（k,k+1:
n）*x（k+1:
n））/A（k,k）;
end
例：
A=[1214;0-42-5;00-5-7.5;000-9];b=[13;2;-35;-18];
x=backsub（A,b）
2、高斯消去法（不选主元素）
functionx=gauss（A,b）
%Input-Aisannxnnonsingularmatrix
%-bisannx1matrix
%Output-xisannx1matrixcontainingthesolutiontoAx=b
%Initializex
[n,n]=size（A）;x=zeros（n,1）;
%Formtheaugmentedmatrix:
Aug=[A|b]
Aug=[Ab];
fork=1:
n-1
%Eliminationprocessforcolumnk
fori=k+1:
n
m=Aug（i,k）/Aug（k,k）;
Aug（i,k:
n+1）=Aug（i,k:
n+1）-m*Aug（k,k:
n+1）;
end
end
%OutputEliminationmatrix
Aug
%BackSubstitution
A=Aug（1:
n,1:
n）;b=Aug（1:
n,n+1）;
x（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（b（k）-A（k,k+1:
n）*x（k+1:
n））/A（k,k）;
end
例：
A=[4-92;2-46;-11-3];b=[5;3;-4];x=gauss（A,b）
3、列主元消去法
functionx=gauss_column（A,b）
%Input-Aisannxnnonsingularmatrix
%-bisannx1matrix
%Output-xisannx1matrixcontainingthesolutiontoAx=b
%InitializexandthetemporarystoragematrixC
[n,n]=size（A）;x=zeros（n,1）;C=zeros（1,n+1）;
%Formtheaugmentedmatrix:
Aug=[A|b]
Aug=[Ab];
fork=1:
n-1
%Partialpivotingforcolumnk
[y,j]=max（abs（Aug（k:
n,k）））;
%Interchangerowkandj
C=Aug（k,:
）;Aug（k,:
）=Aug（j+k-1,:
）;Aug（j+k-1,:
）=C;
ifAug（k,k）==0
'Awassingular,Nouniquesolution'
break
end
%Eliminationprocessforcolumnk
fori=k+1:
n
m=Aug（i,k）/Aug（k,k）;
Aug（i,k:
n+1）=Aug（i,k:
n+1）-m*Aug（k,k:
n+1）;
end
end
%BackSubstitutionon[U|y]usingbacksub
x=backsub（Aug（1:
n,1:
n）,Aug（1:
n,n+1））;
例：
A=[1214;2043;4221;-3132];b=[13;28;20;6];
x=gauss_column（A,b）
4、LU分解法（不选主元素）
functionx=lufact（A,b）
%Input-Aisannxnmatrix
%-bisannx1matrix
%Output-xisannx1matrixcontainingthesolutiontoAx=b
%Initilizex,y,
[n,n]=size（A）;x=zeros（n,1）;y=zeros（n,1）;
fork=1:
n-1
%CalculatemultiplierandplaceinsubdiagonalportionofA
fori=k+1:
n
mult=A（i,k）/A（k,k）;
A（i,k）=mult;
A（i,k+1:
n）=A（i,k+1:
n）-mult*A（k,k+1:
n）;
end
end
%OutputLUfactorizationcompactmatrix
A
%solvefory
y
（1）=b
（1）;
fori=2:
n
y（i）=b（i）-A（i,1:
i-1）*y（1:
i-1）;
end
%Solveforx
x（n）=y（n）/A（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-A（i,i+1:
n）*x（i+1:
n））/A（i,i）;
end
例：
A=[4-92;2-46;-11-3];b=[5;3;-4];x=lufact（A,b）
5、列主元LU分解法（PA＝LU）A是非奇异矩阵
functionx=lu_column（A,b）
%Input-Aisannxnmatrix
%-bisannx1matrix
%Output-xisannx1matrixcontainingthesolutiontoAx=b
%Initilizex,y,thetemporarystoragematrixC,andtherow
%permutationinformationmatrixR
[n,n]=size（A）;
x=zeros（n,1）;y=zeros（n,1）;C=zeros（1,n）;R=1:
n;
fork=1:
n-1
%Findthepivotrowforcolumnk
[maxl,j]=max（abs（A（k:
n,k）））;
%Interchangerowkandj
C=A（k,:
）;A（k,:
）=A（j+k-1,:
）;A（j+k-1,:
）=C;
d=R（k）;R（k）=R（j+k-1）;R（j+k-1）=d;
ifA（k,k）==0
'Aissingular,Nouniquesolution'
break
end
%CalculatemultiplierandplaceinsubdiagonalportionofA
fori=k+1:
n
mult=A（i,k）/A（k,k）;
A（i,k）=mult;
A（i,k+1:
n）=A（i,k+1:
n）-mult*A（k,k+1:
n）;
end
end
%OutputLUfactorizationcompactmatrix
A
%solvefory
y
（1）=b（R
（1））;
fori=2:
n
y（i）=b（R（i））-A（i,1:
i-1）*y（1:
i-1）;
end
%Solveforx
x（n）=y（n）/A（n,n）;
fori=n-1:
-1:
1
x（i）=（y（i）-A（i,i+1:
n）*x（i+1:
n））/A（i,i）;
end
例：
A=[126;48-1;-23-5];b=[1;2;3];
x=lu_column（A,b）
6、Jacobi迭代:
程序可用的充分条件是A具有严格对角优势。
functionx=jacobi（A,b,x0,delta,maxl）
%I
升级会员 