数学数学实验4 Newton迭代法.docx
《数学数学实验4 Newton迭代法.docx》由会员分享,可在线阅读,更多相关《数学数学实验4 Newton迭代法.docx(16页珍藏版)》请在冰豆网上搜索。
![数学数学实验4 Newton迭代法.docx](https://file1.bdocx.com/fileroot1/2022-12/14/2d4a37a9-978b-403e-bd72-11f8556de758/2d4a37a9-978b-403e-bd72-11f8556de7581.gif)
数学数学实验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法寻根的话初值非常不好估计,所以我们要用最稳定的二分法找它们的根。