数学数学实验4 Newton迭代法.docx
- 文档编号:5288850
- 上传时间:2022-12-14
- 格式:DOCX
- 页数:16
- 大小:79.06KB
数学数学实验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; whilek fa=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 iferr r(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法寻根的话初值非常不好估计,所以我们要用最稳定的二分法找它们的根。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学数学实验4 Newton迭代法 数学 实验 Newton 迭代法