Jacobi迭代法Gauss-Seidel迭代法

1、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 - i

2、ndex=1表示迭代收敛到指定要求；%index=0表示迭代失败if nargin <4 it_max=100; endif nargin <3 ep=1e-5; endn=length(A); k=0;x=zeros(n,1); y=zeros(n,1); index=1;while 1for i=1:ny(i)=b(i);for j=1:nif j=iy(i)=y(i)-A(i,j)*x(j);endendif abs(A(i,i)<1e-10 | k=it_maxindex=0; return;endy(i)=y(i)/A(i,i);endif norm(y-x,inf

3、)<epbreak;endx=y; k=k+1;end用Jacobi迭代法求方程组的解。输入：A=4 3 0;3 3 -1;0 -1 4;b=24;30;-24;x, k, index=Jacobi(A, b, 1e-5, 100)输出：x =-2.999811.9987-3.0001k =100index = 02.熟悉Gauss-Seidel迭代法，并编写 Matlab程序function v,sN,vChain=gaussSeidel(A,b,x0,errorBound,maxSp)%Gauss-Seidel迭代法求解线性方程组%A-系数矩阵b-右端向量x0-初始迭代点 error

4、Bound-近似精度maxSp-最大迭代次 数%v-近似解sN-迭代次数vChain-迭代过程的所有值step=0;error=inf;s=size(A);D=zeros(s(1);vChain=zeros(15,3);%最多能记录15次迭代次数k=1;fx0=x0;for i=1:s(1)D(i,i)=A(i,i);end;L=-tril(A,-1);U=-triu(A,1);while error>=errorBound & step<maxSpx0=inv(D)*(L+U)*x0+inv(D)*b;vChain(k,:)=x0'k=k+1;error=norm

5、(x0-fx0);fx0=x0;step=step+1;endv=x0;sN=step;用Gauss-Seidel迭代法求解上题的线性方程组，取 。输入：A=4 3 0;3 3 -1;0 -1 4;b=24;30;-24;x0=0;0;0;v,sN,vChain=gaussSeidel(A,b,x0,0.00001,11)输出：v =0.616911.1962 -4.2056 sN =11 vChain =6.000010.0000-6.0000-1.50002.0000-3.50004.500010.3333-5.5000-1.75003.6667-3.41673.250010.6111-5

6、.0833-1.95835.0556-3.34722.208310.8426-4.7361-2.13196.2130-3.28941.340311.0355-4.4468-2.27667.1775-3.24110.616911.1962-4.2056000000000000例7.6用Gauss-Serdel迭代法求解下列线性方程组。设迭代 初值为0,迭代精度为IO*。在命令中调用函数文件gauseidel.m,命令如下：b=9,7,6，；今xtn=gauseidel(A,b,0,0,01.0e-6)例77分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性 方程组，看是否收敛。命

7、令如下：a=b=9;7;6J;x5n=jacobi(a,b90;0;0)x5n=gauseidel(a9b,0;0;0)数值实验数值实验要求：数值实验报告内容： 要包含题目、算法公式、 完整的程序、正确的数值结果和图形以及相应 的误差分析。在本课程网站上提交数值实验报告的电子文档。一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。1. 对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2. 对这些数据构造Lagrang插值多项式，并画出得到的 Lagrang插值多项式曲线。x 0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.44.75

8、.06.0f(x) 1.3 1.5 1.85 2.1 2.6 2.7 2.42.152.052.1 2.25x 7.0 8.0 9.2 10.5 11.3 11.6 12.012.613.013.3f(x) 2.3 2.25 1.95 1.4 0.9 0.7 0.60.50.40.251.使用三次样条插值函数csape ()求解。解：输入：>>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

9、2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.90.7 0.6 0.5 0.4 0.25;>> pp=csape(x,y,'variational',0 0)得到：pp =form: 'pp'breaks: 0.9000 1.3000 1.9000 2.1000 2.6000 3 3.9000 4.4000 4.7000 5 6 7 8 9.2000 10.500011.3000 11.6000 12 12.6000 13 13.3000coefs: 20x4 doublepieces: 20o

10、rder: 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

11、.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=i

12、nterp1(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

13、,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*xA20+.42770e-8*xA19-.27708e-6*xA18+.11098e-4*xA17-.30784e-3*xA16+.62785 e-2*xA15-.97558e-1*xA14+1.1810*xA13-11.297*xA12+86.085*xA11-524.68*xA10+2558

14、.0*xA9-9 942.3*xA8+30586.*xA7-73622.*xA6+.13627e6*xA5-.18907e6*xA4+.18913e6*xA3-.12803e6*xA2+52170.*x-9593.4再输入：>> y1=polyval(a,x);plot(x,y1);xlabel('x');ylabel('y');得到图形：结果分析：由以上两图形可以看出三次样条插值的图形较lagrange插值的图形要光滑的多。因为样条函数插值不仅具有较好的收敛性和稳定性，而且其光滑性也较高，因此样条函数成为了重要的插值工具，其中应用较多的是三次样条插

15、值。二、对于问题将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构成的向

16、量， y为Yn构成的向量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=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-xA2+1','x','y');t,u=Euler(ydot,0,0.5,0.025,20)得到数值结果：t =Columns 1 through 1700.02500.05000.07500.10000.12500.150

17、00.17500.20000.22500.25000.27500.30000.32500.35000.37500.4000Columns 18 through 210.42500.45000.47500.5000 u =Columns 1 through 170.50000.53750.57590.61530.65550.69660.73870.78160.82530.87000.91550.96181.00891.05691.10571.15531.2056Columns 18 through 211.25681.30871.36131.4147即采用Euler法得到：u(0.1)=0.65

18、55,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;for n

19、=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-xA2+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.5768

20、0.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构成的向

21、量x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=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-xA2+1'

22、;,'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结果比较：t Euler法改进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.42

23、501.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

24、|< §则停止计算。(4) 若k=N，则停止计算；否则置 k=k+1，转(2)。Matlab 程序：function x_star,index,it= Newton(fun,x,ep,it_max)%求解非线性方程的Newton法% fun (x)为需要求根的函数，第一个分量是函数值，第二个分量是导数值% fun=inline('xA3-x-1,3*xA2-1')当 f(x)=xA3-x-1;%x为初始点% ep为精度，缺省值为 1e-5% it_max为最大迭代次数，缺省值 100% x_star为当迭代成功时输出方程的根，失败时输出最后的迭代值% index

25、为指标变量，index=1表示迭代成功index=0表示失败% it为迭代次数if nargin<4it_max=100;endif nargin<3ep=1e-5;endindex=0;k=1;while k<=it_maxx1=x; f=feval(fun,x);if abs(f(2)<ep break; endx=x-f(1)/f(2);if abs(x-x1)<ep index=1;break;endk=k+1;endx_star=x;it=k;解：(1)由于 f(x)=arctan(x) , f ' (x)= 1/1+x2!X初始值 x0=1.4

26、5,输入>> fun=inline('atan(x),1/(1+xA2)');>> x_star,index,it=Newton(fun,1.45)得到数值结果：x_star =1.6281e+004index = 0it =7取初始值x0=0.5，输入>> fun=inline('atan(x),1/(1+xA2)');>> x_star,index,it=Newton(fun,0.5)得到数值结果：x_star =0index = 1it =4输入 x=-3:0.05:3; y=atan(x);plot(x,y

27、);xlabel('x');ylabel('y');得到图形：(2)由于 f(x)=x3-x- 3=0, f '=(班2-1,取初始值 x0=0.0,输入>> fun=inline('xA3-x-3,3*xA2-1');>> x_star,index,it=Newton(fun,0.0)得到数值结果：x_star = -0.0074index = 0it =101取初始值x0=2.0，输入>> fun=inline('xA3-x-3,3*xA2-1');>> x_star,i

28、ndex,it=Newton(fun,2.0)得到数值结果：x_star = 1.6717index = 1it =4输入 x=0:0.05:3; y=x.A3-x-3;plot(x,y);xlabel('x');ylabel('y');得到图形：结果分析：从以上结果可以看出初值的选取与Newton法的收敛很有关系。初值选的不好，Nexton法将不收敛，它的收敛性是在跟 a附近讨论的。选取初值时一定要十分小心。四、用Jacobi迭代和SOR法分别求解线性方程组（教科书第86页算例2）,并验证、输出SOR法的松弛因子 w和对应的迭代次数。1. Jacobi迭代法J

29、acobi迭代法的算法为Matlab 程序：functionx,k,index=Jacobi(A,b,ep,it_max)%求解线性方程组的Jacobi迭代法% A为系数矩阵%b为方程组右端项% ep为精度要求，缺省值 1e-5% it_max为最大迭代次数，缺省值 100%x为方程组的解% k为迭代次数% index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。if nargin<4 it_max=100;endif nargin<3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=

30、1;while 1for i=1:ny(i)=b(i);for j=1:nif j=iy(i)=y(i)-A(i,j)*x(j);endendif abs(A(i,i)<1e-10|k=it_maxindex=0;return;endy(i)=y(i)/A(i,i);endif norm(y-x,inf)<epbreak;endx=y;k=k+1;endfunction A,b=shuru%自己定义输入矩阵A和向量b的函数n=15;for i=1:nif i=1|i=15z(i)=3;else z(i)=2;endfor j=1:nif j=i A(i,j)=4;elseif j=

31、i+1|j=i-1 A(i,j)=-1;else A(i,j)=0;endendendb=z'解：输入：>> A,b=shuru;ep=1e-6;>> x,k,index=Jacobi(A,b,ep)得到数值结果：x =(1.0000,1.0000,- ,1.0000)Tk = 19index =12. so*SOR迭代法的算法为:Matlab 程序：function x,k,index=SOR1(A,b,ep,w,it_max)%求解线性方程组的 SOR迭代法% A为系数矩阵%b为方程组右端项% ep为精度要求，缺省值 1e-5% w为超松弛因子，缺省值为1;

32、% it_max为最大迭代次数，缺省值 100%x为方程组的解% k为迭代次数% index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。if nargin<5 it_max=100;endif nargin<4 w=1;endif nargin<3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1y=x;for i=1:nz=b(i);for j=1:nif j=iz=z-A(i,j)*x(j);endendif abs(A(i,i)<1e-10|k=i

33、t_max index=0;return;endz=z/A(i,i);x(i)=(1-w)*x(i)+w*z;endif norm(y-x,inf)<epbreak;endk=k+1;end解：取w=1.1，输入：>> 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.8 190.9 161.0 131.1 1

34、11.2 141.3 171.4 231.5 27结果分析:从以上结果可以看出在求解相同问题时，sor法的迭代次数要比Jacobi迭代法少很多，计算量小许多。此外可以看出松弛因子w的选取对迭代次数的影响也十分大。在实际计算时，最优松弛因子很难事先确定，一般可用试算法取近似最优值。是n维非奇异阵是n维向量 是初值delta是误差界X为所求的方程组 AX=B的近似解k=1:max1j=1:N(err<delta)A=4,1,-1;1,-5,-1;2,-1,-6>>P=0;0;0>>X=jacobi(A,b,P,10A(-4),20)=nargin是用来判断输入变量个数

35、的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值，如下面的函数。例子，函数test 1的功能是输出a和b的和。如果只输入一个变量，则认为另一个变量为0,如果两个变量都没有输入，则默认两者均为0 。nargin=0a=0;b=0;elseif nargin=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+