jacobi迭代法 gauss-seidel迭代法

《jacobi迭代法 gauss-seidel迭代法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

Matlab线性方程组的迭代解法（Jacobi迭代法Gauss-Seidel迭代法）实验报告2008年11月09日星期日12:491.熟悉Jacobi迭代法，并编写Matlab程序matlab程序按照算法(Jacobi迭代法)编写Matlab程序(Jacobi.m)function[x,k,index]=Jacobi(A,b,ep,it_max)%求解线性方程组的Jacobi迭代法，其中%A---方程组的系数矩阵%b---方程组的右端项%ep---精度要求。省缺为1e-5%it_max---最大迭代次数，省缺为100%x---方程组的解%k---迭代次数%index---index=1表示迭代收敛到指定要求；%index=0表示迭代失败ifnargin<4it_max=100;endifnargin<3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1fori=1:ny(i)=b(i);forj=1:nifj~=iy(i)=y(i)-A(i,j)*x(j);endendifabs(A(i,i))<1e-10|k==it_maxindex=0;return;endy(i)=y(i)/A(i,i);endifnorm(y-x,inf)=errorBound&step>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.612.012.613.013.3];>>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];>>pp=csape(x,y,'variational',[00])得到：pp=form:'pp'breaks:[0.90001.30001.90002.10002.600033.90004.40004.700056789.200010.500011.300011.60001212.60001313.3000]coefs:[20x4double]pieces:20order:4dim:1再输入：>>pp.coefs得到：ans=-0.247600.53961.30000.9469-0.29720.42081.5000-2.95641.40731.08681.8500-0.4466-0.36661.29492.10000.4451-1.03650.59342.60000.1742-0.5025-0.02222.70000.0781-0.0322-0.50342.40001.31420.0849-0.47712.1500-1.58121.2676-0.07132.05000.0431-0.15550.26232.1000-0.0047-0.02610.08082.2500 -0.0244-0.04010.01462.30000.0175-0.1135-0.13902.2500-0.0127-0.0506-0.33581.9500-0.0203-0.1002-0.53181.40001.2134-0.1490-0.73120.9000-0.83930.9431-0.49290.70000.0364-0.0640-0.14130.6000-0.44800.0014-0.17890.50000.5957-0.5361-0.39280.4000因此，三次样条函数S（x）为最后输入：>>xx=0:0.01:14;yy=interp1(x,y,xx,'csape');>>plot(xx,yy)；xlabel('x');ylabel('y');得到图形：2.Lagrange插值算法：1)输入N个节点数组；2)定义初始值和；3)利用公式求N次插值基函数；4)利用多项式求插值函数;解：输入：>>x=[0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3];y=[1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25];>>a=polyfit(x,y,length(x)-1);>>p=vpa(poly2sym(a),5)得到数值结果：p=.30732e-10*x^20+.42770e-8*x^19-.27708e-6*x^18+.11098e-4*x^17-.30784e-3*x^16+.62785e-2*x^15-.97558e-1*x^14+1.1810*x^13-11.297*x^12+86.085*x^11-524.68*x^10+2558.0*x^9-9942.3*x^8+30586.*x^7-73622.*x^6+.13627e6*x^5-.18907e6*x^4+.18913e6*x^3-.12803e6*x^2+52170.*x-9593.4再输入：>>y1=polyval(a,x);plot(x,y1);xlabel('x');ylabel('y'); 得到图形：结果分析：由以上两图形可以看出三次样条插值的图形较lagrange插值的图形要光滑的多。因为样条函数插值不仅具有较好的收敛性和稳定性，而且其光滑性也较高，因此样条函数成为了重要的插值工具，其中应用较多的是三次样条插值。二、对于问题将h=0.025的Euler法，h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0.1，0.2，0.3，0.4和0.5处进行比较。精确解为：。1.Euler法在x,y平面上微分方程的解在曲线上一点(x,y)的切线斜率等于函数的值。该曲线的顶点设为p,再推进到p(),显然两个顶点p,p的坐标有以下关系Matlab程序：function[x,y]=Euler(ydot,x0,y0,h,N)%ydot为一阶微分方程的函数%x0，y0为初始条件%h为区间步长%N为区间个数%x为Xn构成的向量，y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(ydot,x(n),y(n));end解：取h=0.025，N=20，输入：>>ydot=inline('y-x^2+1','x','y');[t,u]=Euler(ydot,0,0.5,0.025,20)得到数值结果：t=Columns1through1700.02500.05000.07500.10000.12500.15000.17500.2000 0.22500.25000.27500.30000.32500.35000.37500.4000Columns18through210.42500.45000.47500.5000u=Columns1through170.50000.53750.57590.61530.65550.69660.73870.78160.82530.87000.91550.96181.00891.05691.10571.15531.2056Columns18through211.25681.30871.36131.4147即采用Euler法得到：u(0.1)=0.6555,u(0.2)=0.8253,u(0.3)=1.0089,u(0.4)=1.2056,u(0.5)=1.41472.改进Euler方法改进Euler公式.y=yn+h/2(f()+f(xn+1,+h*f()))即迭代公式为:Matlab程序：function[x,y]=Euler_ad(ydot,x0,y0,h,N)%改进Euler公式%ydot为一阶微分方程的函数%x0，y0为初始条件%h为区间步长%N为区间个数%x为Xn构成的向量，y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;ybar=y(n)+h*feval(ydot,x(n),y(n));y(n+1)=y(n)+h/2*(feval(ydot,x(n),y(n))+feval(ydot,x(n+1),ybar));end解：取h=0.05，N=10，输入： >>ydot=inline('y-x^2+1','x','y');[t,u]=Euler_ad(ydot,0,0.5,0.05,10)得到数值结果：t=00.05000.10000.15000.20000.25000.30000.35000.40000.45000.5000u=0.50000.57680.65730.74140.82910.92021.01471.11261.21361.31781.4250即采用改进的Euler法得到：u(0.1)=0.6573,u(0.2)=0.8291,u(0.3)=1.0147,u(0.4)=1.2136,u(0.5)=1.42503.四阶Runge_Kutta法求解公式为：Matlab程序：function[x,y]=Runge_Kutta4(ydot,x0,y0,h,N)%标准四阶Runge_Kutta方法%ydot为一阶微分方程的函数%x0，y0为初始条件%h为区间步长%N为区间个数%x为Xn构成的向量，y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;k1=h*feval(ydot,x(n),y(n));k2=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k1);k3=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k2);k4=h*feval(ydot,x(n)+h,y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end 解：取h=0.1，N=5，输入：>>ydot=inline('y-x^2+1','x','y');[t,u]=Runge_Kutta4(ydot,0,0.5,0.1,5)得到数值结果：t=00.10000.20000.30000.40000.5000u=0.50000.65740.82931.01511.21411.4256结果比较：tEuler法改进Euler法四阶Runge_Kutta精确解0.10.65550.65730.65740.65740.20.82530.82910.82930.82930.31.00891.01471.01511.01510.41.20561.21361.21411.21410.51.41471.42501.42561.4256由以上结果可以看出改进的Euler法较Euler法计算精度有所提高，但还不是十分精确。四阶Runge_Kutta法具有非常高的精度，事实上，在求解微分方程初值问题，四阶Runge_Kutta法是单步长中最优秀的方法，通常都用该方法求解实际问题。三、用Newton迭代法求方程的根时，分别取初始值，；用Newton迭代法求方程时，分别取初始值，；算法：（1）取初始点x0最大迭代次数N和精度要求ε，k=0.（2）如果f’（xk）=0，则停止计算；否则计算Xk+1=xk-f（xk）/f’（xk）（3）若|xk+1-xk|<ε，则停止计算。（4）若k=N，则停止计算；否则置k=k+1，转（2）。Matlab程序：function[x_star,index,it]=Newton(fun,x,ep,it_max)%求解非线性方程的Newton法%fun（x)为需要求根的函数，第一个分量是函数值，第二个分量是导数值%fun=inline('[x^3-x-1,3*x^2-1]')当f(x)=x^3-x-1;%x为初始点%ep为精度，缺省值为1e-5%it_max为最大迭代次数，缺省值100%x_star为当迭代成功时输出方程的根，失败时输出最后的迭代值%index为指标变量，index=1表示迭代成功index=0表示失败%it为迭代次数 ifnargin<4it_max=100;endifnargin<3ep=1e-5;endindex=0;k=1;whilek<=it_maxx1=x;f=feval(fun,x);ifabs(f(2))>fun=inline('[atan(x),1/(1+x^2)]');>>[x_star,index,it]=Newton(fun,1.45)得到数值结果：x_star=1.6281e+004index=0it=7取初始值x0=0.5,输入>>fun=inline('[atan(x),1/(1+x^2)]');>>[x_star,index,it]=Newton(fun,0.5)得到数值结果：x_star=0index=1it=4输入x=-3:0.05:3;y=atan(x);plot(x,y);xlabel('x');ylabel('y');得到图形：(2)由于f(x)=x3-x-3=0,f’(x)=3x2-1,取初始值x0=0.0,输入>>fun=inline('[x^3-x-3,3*x^2-1]');>>[x_star,index,it]=Newton(fun,0.0)得到数值结果：x_star=-0.0074index=0it=101取初始值x0=2.0,输入>>fun=inline('[x^3-x-3,3*x^2-1]');>>[x_star,index,it]=Newton(fun,2.0) 得到数值结果：x_star=1.6717index=1it=4输入x=0:0.05:3;y=x.^3-x-3;plot(x,y);xlabel('x');ylabel('y');得到图形：结果分析：从以上结果可以看出初值的选取与Newton法的收敛很有关系。初值选的不好，Nexton法将不收敛，它的收敛性是在跟a附近讨论的。选取初值时一定要十分小心。四、用Jacobi迭代和SOR法分别求解线性方程组（教科书第86页算例2），并验证、输出SOR法的松弛因子w和对应的迭代次数。1.Jacobi迭代法Jacobi迭代法的算法为:Matlab程序：function[x,k,index]=Jacobi(A,b,ep,it_max)%求解线性方程组的Jacobi迭代法%A为系数矩阵%b为方程组右端项%ep为精度要求，缺省值1e-5%it_max为最大迭代次数，缺省值100%x为方程组的解%k为迭代次数%index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。ifnargin<4it_max=100;endifnargin<3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1fori=1:ny(i)=b(i);forj=1:nifj~=i y(i)=y(i)-A(i,j)*x(j);endendifabs(A(i,i))<1e-10|k==it_maxindex=0;return;endy(i)=y(i)/A(i,i);endifnorm(y-x,inf)>[A,b]=shuru;ep=1e-6;>>[x,k,index]=Jacobi(A,b,ep)得到数值结果：x=(1.0000,1.0000,…,1.0000)Tk=19index=1 2.sor法SOR迭代法的算法为：Matlab程序：function[x,k,index]=SOR1(A,b,ep,w,it_max)%求解线性方程组的SOR迭代法%A为系数矩阵%b为方程组右端项%ep为精度要求，缺省值1e-5%w为超松弛因子，缺省值为1；%it_max为最大迭代次数，缺省值100%x为方程组的解%k为迭代次数%index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。ifnargin<5it_max=100;endifnargin<4w=1;endifnargin<3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while1y=x;fori=1:nz=b(i);forj=1:nifj~=iz=z-A(i,j)*x(j);endendifabs(A(i,i))<1e-10|k==it_maxindex=0;return;endz=z/A(i,i);x(i)=(1-w)*x(i)+w*z;endifnorm(y-x,inf)>[A,b]=shuru;ep=1e-6;w=1.1; >>[x,k,index]=sor(A,b,ep,w)得到数值结果：x=(1.0000,1.0000,…,1.0000)Tk=11index=1依次取w=0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5得到下表：松弛因子w迭代次数0.8190.9161.0131.1111.2141.3171.4231.527结果分析：从以上结果可以看出在求解相同问题时，sor法的迭代次数要比Jacobi迭代法少很多，计算量小许多。此外可以看出松弛因子w的选取对迭代次数的影响也十分大。在实际计算时，最优松弛因子很难事先确定，一般可用试算法取近似最优值。functionX=jacobi(A,b,P,delta,max1)%A是n维非奇异阵%B是n维向量%P是初值%delta是误差界%X为所求的方程组AX=B的近似解N=length(b);fork=1:max1forj=1:NX(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));P=X';if(err>A=[4,1,-1;1,-5,-1;2,-1,-6]>>b=[13;-8;-2]>>P=[0;0;0]>>X=jacobi(A,b,P,10^(-4),20)k=9err=2.5713e-005X=3.00002.00001.0000nargin是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值，如下面的函数。例子，函数test1的功能是输出a和b的和。如果只输入一个变量，则认为另一个变量为0，如果两个变量都没有输入，则默认两者均为0。functiony=test1(a,b)ifnargin==0a=0;b=0;elseifnargin==1b=0;endy=a+b;s （2）迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。i．Jacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。Jacobi迭代法的MATLAB函数文件Jacobi.m如下：function[y,n]=jacobi(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D(L+U);f=Db;y=B*x0+f;n=1;%迭代次数whilenorm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end例用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在命令中调用函数文件Jacobi.m，命令如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)ii．Gauss-Serdel迭代法在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel 迭代用新分量代替旧分量，精度会高些。Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function[y,n]=gauseidel(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)U;f=(D-L)b;y=G*x0+f;n=1;%迭代次数whilenorm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end例用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在命令中调用函数文件gauseidel.m，命令如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)例分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。命令如下：a=[1,2,-2;1,1,1;2,2,1];b=[9;7;6];[x,n]=jacobi(a,b,[0;0;0])[x,n]=gauseidel(a,b,[0;0;0])