数值方法.docx
- 文档编号:4325005
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:42
- 大小:220.21KB
数值方法.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(err end 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-a end 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(err end 例: 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(err end 例: 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(err break 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 方法