数学数学实验4 Newton迭代法.docx
《数学数学实验4 Newton迭代法.docx》由会员分享，可在线阅读，更多相关《数学数学实验4 Newton迭代法.docx（16页珍藏版）》请在冰豆网上搜索。
数学数学实验4Newton迭代法
数学实验题目4Newton迭代法
摘要
为初始猜测，则由递推关系
产生逼近解
的迭代序列
，这个递推公式就是Newton法。
当
距
较近时，
很快收敛于
。
但当
选择不当时，会导致
发散。
故我们事先规定迭代的最多次数。
若超过这个次数，还不收敛，则停止迭代另选初值。
前言
利用牛顿迭代法求
的根
程序设计流程
问题1
（1）程序运行如下：
r=NewtSolveOne（'fun1_1',pi/4,1e-6,1e-4,10）
r=0.7391
（2）程序运行如下：
r=NewtSolveOne（'fun1_2',0.6,1e-6,1e-4,10）
r=0.5885
问题2
（1）程序运行如下：
r=NewtSolveOne（'fun2_1',0.5,1e-6,1e-4,10）
r=0.5671
（2）程序运行如下：
r=NewtSolveOne（'fun2_2',0.5,1e-6,1e-4,20）
r=0.5669
问题3
（1）程序运行如下：
①p=LegendreIter
（2）
p=1.00000-0.3333
p=LegendreIter（3）
p=1.00000-0.60000
p=LegendreIter（4）
p=1.00000-0.857100.0857
p=LegendreIter（5）
p=1.00000-1.111100.23810
②p=LegendreIter（6）
p=1.00000-1.363600.45450-0.0216
r=roots（p）'
r=-0.932469514203150-0.6612093864662650.932469514203153
0.661209386466264-0.2386191860831970.238619186083197
用二分法求根为：
r=BinSolve（'LegendreP6',-1,1,1e-6）
r=-0.932470204878826-0.661212531887755-0.238620057397959
0.2386001275510200.6611926020408160.932467713647959
（2）程序运行如下：
①p=ChebyshevIter
（2）
p=1.00000-0.5000
p=ChebyshevIter（3）
p=1.00000-0.75000
p=ChebyshevIter（4）
p=1.00000-1.000000.1250
p=ChebyshevIter（5）
p=1.00000-1.250000.31250
②p=ChebyshevIter（6）
p=1.00000-1.500000.56250-0.0313
r=roots（p）'
r=-0.965925826289067-0.7071067811865480.965925826289068
0.707106781186547-0.2588190451025210.258819045102521
用二分法求根为：
r=BinSolve（'ChebyshevT6',-1,1,1e-6）
r=-0.965929926658163-0.707110969387755-0.258828922193878
0.2588189572704080.7071059869260200.965924944196429
与下列代码结果基本一致，只是元素顺序稍有不同：
j=0:
5;
x=cos（（2*j+1）*pi/2/（5+1））
x=0.9659258262890680.7071067811865480.258819045102521
-0.258819045102521-0.707106781186547-0.965925826289068
（3）程序运行如下：
①p=LaguerreIter
（2）
p=1-42
p=LaguerreIter（3）
p=1-918-6
p=LaguerreIter（4）
p=1-1672-9624
p=LaguerreIter（5）
p=1.0000-25.0000200.0000-600.0000600.0000-120.000
②p=LaguerreIter（5）
p=1.0000-25.0000200.0000-600.0000600.0000-120.000
r=roots（p）'
r=12.6408008442757327.0858100058588913.596425771040711
1.4134030591065200.263560319718141
用二分法求根为：
r=BinSolve（'LaguerreL5',0,13,1e-6）
r=0.2635603145677221.4134030561057893.596425765631150
.0858********
（4）程序运行如下：
①p=HermiteIter
（2）
p=1.00000-0.5000
p=HermiteIter（3）
p=1.00000-1.50000
p=HermiteIter（4）
p=1.00000-3.000000.7500
p=HermiteIter（5）
p=1.00000-5.000003.75000
②p=HermiteIter（6）
p=1.00000-7.5000011.25000-1.8750
r=roots（p）'
r=-2.3506049736744872.350604973674488-1.335849074013696
1.335849074013698-0.4360774119276170.436077411927616
用二分法求根为：
r=BinSolve（'HermiteH6',-3,3,1e-6）
r=-2.350604981792216-1.335849*********-0.436077818578603
0.4360773514728161.3358489834532442.350604952598104
所用到的函数
functionr=NewtSolveOne（fun,x0,ftol,dftol,maxit）
%NewtSolveOne用Newton法解方程f（x）=0在x0附近的一个根
%
%Synopsis:
r=NewtSolveOne（fun,x0）
%r=NewtSolveOne（fun,x0,ftol,dftol）
%
%Input:
fun=（string）需要求根的函数及其导数
%x0=猜测根，Newton法迭代初始值
%ftol=（optional）误差，默认为5e-9
%dftol=（optional）导数容忍最小值，小于它表明Newton法失败，默认为5e-9
%maxit=（optional）迭代次数，默认为25
%
%Output:
r=在寻根区间内的根或奇点
ifnargin<3
ftol=5e-9;
end
ifnargin<4
dftol=5e-9;
end
ifnargin<5
maxit=25;
end
x=x0;%设置初始迭代位置为x0
k=0;%初始化迭代次数为0
whilek<=maxit
k=k+1;
[f,dfdx]=feval（fun,x）;%fun返回f（x）和f'（x）的值
ifabs（dfdx）r=[];
warning（'dfdxistoosmall!
'）;
return;
end
dx=f/dfdx;%x（n+1）=x（n）-f（x（n））/f'（x（n）），这里设dx=f（x（n））/f'（x（n））
x=x-dx;
ifabs（f）r=x;
return;
end
end
r=[];%如果牛顿法未收敛，返回空值
functionp=LegendreIter（n）
%LegendreIter用递推的方法计算n次勒让德多项式的系数向量Pn+2（x）=（2*i+3）/（i+2）*x*Pn+1（x）-（i+1）/（i+2）*Pn（x）
%
%Synopsis:
p=LegendreIter（n）
%
%Input:
n=勒让德多项式的次数
%
%Output:
p=n次勒让德多项式的系数向量
ifround（n）~=n|n<0
error（'n必须是一个非负整数'）;
end
ifn==0%P0（x）=1
p=1;
return;
elseifn==1%P1（x）=x
p=[10];
return;
end
pBk=1;%初始化三项递推公式后项为P0
pMid=[10];%初始化三项递推公式中项为P1
fori=0:
n-2
pMidCal=zeros（1,i+3）;%构造用于计算的x*Pn+1
pMidCal（1:
i+2）=pMid;
pBkCal=zeros（1,i+3）;%构造用于计算的Pn
pBkCal（3:
i+3）=pBk;
pFwd=（2*i+3）/（i+2）*pMidCal-（i+1）/（i+2）*pBkCal;%勒让德多项式三项递推公式Pn+2（x）=（2*i+3）/（i+2）*x*Pn+1（x）-（i+1）/（i+2）*Pn（x）
pBk=pMid;%把中项变为后项进行下次迭代
pMid=pFwd;%把前项变为中项进行下次迭代
end
p=pFwd/pFwd
（1）;%把勒让德多项式最高次项系数归一化
functionp=ChebyshevIter（n）
%ChebyshevIter用递推的方法计算n次勒让德-切比雪夫多项式的系数向量Tn+2（x）=2*x*Tn+1（x）-Tn（x）
%
%Synopsis:
p=ChebyshevIter（n）
%
%Input:
n=勒让德-切比雪夫多项式的次数
%
%Output:
p=n次勒让德-切比雪夫多项式的系数向量
ifround（n）~=n|n<0
error（'n必须是一个非负整数'）;
end
ifn==0%T0（x）=1
p=1;
return;
elseifn==1%T1（x）=x
p=[10];
return;
end
pBk=1;%初始化三项递推公式后项为T0
pMid=[10];%初始化三项递推公式中项为T1
fori=0:
n-2
pMidCal=zeros（1,i+3）;%构造用于计算的x*Tn+1
pMidCal（1:
i+2）=pMid;
pBkCal=zeros（1,i+3）;%构造用于计算的Pn
pBkCal（3:
i+3）=pBk;
pFwd=2*pMidCal-pBkCal;%勒让德-切比雪夫多项式三项递推公式Tn+2（x）=2*x*Tn+1（x）-Tn（x）
pBk=pMid;%把中项变为后项进行下次迭代
pMid=pFwd;%把前项变为中项进行下次迭代
end
p=pFwd/pFwd
（1）;%把勒让德-切比雪夫多项式最高次项系数归一化
functionp=LaguerreIter（n）
%LaguerreIter用递推的方法计算n次拉盖尔多项式的系数向量Ln+2（x）=（2*n+3-x）*Ln+1（x）-（n+1）*Ln（x）
%
%Synopsis:
p=LaguerreIter（n）
%
%Input:
n=拉盖尔多项式的次数
%
%Output:
p=n次拉盖尔多项式的系数向量
ifround（n）~=n|n<0
error（'n必须是一个非负整数'）;
end
ifn==0%L0（x）=1
p=1;
return;
elseifn==1%L1（x）=-x+1
p=[-11];
return;
end
pBk=1;%初始化三项递推公式后项为L0
pMid=[-11];%初始化三项递推公式中项为L1
fori=0:
n-2
pMidCal1=zeros（1,i+3）;%构造用于计算的x*Ln+1（x）
pMidCal1（1:
i+2）=pMid;
pMidCal2=zeros（1,i+3）;%构造用于计算的Ln+1（x）
pMidCal2（2:
i+3）=pMid;
pBkCal=zeros（1,i+3）;%构造用于计算的Ln（x）
pBkCal（3:
i+3）=pBk;
pFwd=（（2*i+3）*pMidCal2-pMidCal1-（i+1）*pBkCal）/（i+2）;%拉盖尔多项式三项递推公式Ln+2（x）=（2*n+3-x）*Ln+1（x）-（n+1）^2*Ln（x）
pBk=pMid;%把中项变为后项进行下次迭代
pMid=pFwd;%把前项变为中项进行下次迭代
end
p=pFwd/pFwd
（1）;%把拉盖尔多项式最高次项系数归一化
functionp=HermiteIter（n）
%HermiteIter用递推的方法计算n次埃尔米特多项式的系数向量Hn+2（x）=2*x*Hn+1（x）-2*（n+1）*Hn（x）
%
%Synopsis:
p=HermiteIter（n）
%
%Input:
n=埃尔米特多项式的次数
%
%Output:
p=n次埃尔米特多项式的系数向量
ifround（n）~=n|n<0
error（'n必须是一个非负整数'）;
end
ifn==0%H0（x）=1
p=1;
return;
elseifn==1%H1（x）=2*x
p=[20];
return;
end
pBk=1;%初始化三项递推公式后项为L0
pMid=[20];%初始化三项递推公式中项为L1
fori=0:
n-2
pMidCal=zeros（1,i+3）;%构造用于计算的x*Hn+1（x）
pMidCal（1:
i+2）=pMid;
pBkCal=zeros（1,i+3）;%构造用于计算的Hn（x）
pBkCal（3:
i+3）=pBk;
pFwd=2*pMidCal-2*（i+1）*pBkCal;%埃尔米特多项式三项递推公式Hn+2（x）=2*x*Hn+1（x）-2*（n+1）*Hn（x）
pBk=pMid;%把中项变为后项进行下次迭代
pMid=pFwd;%把前项变为中项进行下次迭代
end
p=pFwd/pFwd
（1）;%把拉盖尔多项式最高次项系数归一化
functionr=BinSolve（fun,a,b,tol）
%BinSolve用二分法解方程f（x）=0在区间[a,b]的根
%
%Synopsis:
r=BinSolve（fun,a,b）
%r=BinSolve（fun,a,b,tol）
%
%Input:
fun=（string）需要求根的函数
%a,b=寻根区间上下限
%tol=（optional）误差，默认为5e-9
%
%Output:
r=在寻根区间内的根
ifnargin<4
tol=5e-9;
end
Xb=RootBracket（fun,a,b）;%粗略寻找含根区间
[m,n]=size（Xb）;
r=[];
nr=1;%初始化找到的根的个数为1
maxit=50;%最大二分迭代次数为50
fori=1:
m
a=Xb（i,1）;%初始化第i个寻根区间下限
b=Xb（i,2）;%初始化第i个寻根区间上限
err=1;%初始化误差
k=0;
whilekfa=feval（fun,a）;%计算下限函数值
fb=feval（fun,b）;%计算上限函数值
m=（a+b）/2;
fm=feval（fun,m）;
err=abs（fm）;
ifsign（fm）==sign（fb）%若中点处与右端点函数值同号，右端点赋值为中点
b=m;
else%若中点处与左端点函数值同号或为0，左端点赋值为中点
a=m;
end
iferrr（nr）=a;%一般奇点不符合该条件，这样可以去除奇点
nr=nr+1;%找到根的个数递增
k=maxit;%改变k值跳出循环
end
k=k+1;%二分迭代次数递增
end
end
functionX=powerX（x,a,b）
%powerX对给定向量（x1,x2,...,xn）返回增幂矩阵（x1^a,x2^a,...,xn^a;x1^a+1,x2^a+1,...,xn^a+1;...;x1^b,x2^b,...,xn^b;）
%
%Synopsis:
X=powerX（x,a,b）
%
%Input:
x=需要返回增幂矩阵的向量
%a,b=寻根区间上下限
%
%Output:
X=增幂矩阵（x1^a,x2^a,...,xn^a;x1^a+1,x2^a+1,...,xn^a+1;...;x1^b,x2^b,...,xn^b;）
ifround（a）~=a|round（b）~=b
error（'a,bmustbeintegers'）;
elseifa>=b
error（'amustbesmallerthanb!
'）;
end
x=x（:
）';
row=b-a+1;
col=length（x）;
X=zeros（row,col）;
fori=b:
-1:
a
X（b-i+1,:
）=x.^i;
End
function[f,dfdx]=fun1_1（x）
f=cos（x）-x;
dfdx=-sin（x）-1;
function[f,dfdx]=fun1_2（x）
f=exp（-x）-sin（x）;
dfdx=-exp（-x）-cos（x）;
function[f,dfdx]=fun2_1（x）
f=x-exp（-x）;
dfdx=1+exp（-x）;
function[f,dfdx]=fun2_2（x）
f=x.^2-2*x*exp（-x）+exp（-2*x）;
dfdx=2*x-2*exp（-x）+2*x*exp（-x）-2*exp（-2*x）;
functiony=LegendreP6（x）
p=LegendreIter（6）;
X=powerX（x,0,6）;
y=p*X;
functiony=ChebyshevT6（x）
p=ChebyshevIter（6）;
X=powerX（x,0,6）;
y=p*X;
functiony=LaguerreL5（x）
p=LaguerreIter（5）;
X=powerX（x,0,5）;
y=p*X;
functiony=HermiteH6（x）
p=HermiteIter（6）;
X=powerX（x,0,6）;
y=p*X;
思考题
（1）由于Newton法具有局部收敛性，所以在实际问题中，当实际问题本身能提供接近于根的初始近似值时，就可保证迭代序列收敛，但当初值难以确定时，迭代序列就不一定收敛。
实际计算时应先用比较稳定的算法，如二分法，计算根的近似值，再将该近似值作为牛顿法的初值，以保证迭代序列的收敛性。
（2）实验2中两个方程根其实相同，只是第二个方程为重根，通过比较迭代次数，第一个方程迭代了3次得出结果，第二个方程迭代了8次得出结果，且第二个方程的结果不如第一个准确，这是由于第二个方程在根处导数为0，在根的领域内导数很小使Newton法收敛速度变慢，精度变低。
（3）我们来看下这几个多项式的图形：
LegendreP6ChebyshevT6
LaguerreL5HermiteH6
我们发现，这些多项式在比较小的区间内有多个根，这就致使其导数也会有多个根，因此如果用Newton法寻根的话初值非常不好估计，所以我们要用最稳定的二分法找它们的根。
升级会员 