Jacobi迭代法_Gauss-Seidel迭代法.docVIP

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表示迭代失败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 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end if abs(A(i,i)1e-10 | k=it_max index=0; return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf)=errorBound & stepx=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; 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.5000 11.3000 11.6000 12 12.6000 13 13.3000 coefs: 20x4 double pieces: 20 order: 4 dim: 1 再输入： pp.coefs 得到： ans = -0.2476 0 0.5396 1.3000 0.9469 -0.2972 0.4208 1.5000 -2.9564 1.4073 1.0868 1.8500 -0.4466 -0.3666 1.2949 2.1000 0.4451 -1.0365 0.5934 2.6000 0.1742 -0.5025 -0.0222 2.7000 0.0781 -0.0322 -0.5034 2.4000 1.3142 0.0849 -0.4771 2.1500 -1.5812 1.2676 -0.0713 2.0500 0.0431 -0.1555 0.2623 2.1000 -0.0047 -0.0261 0.0808 2.2500 -0.0244 -0.0401 0.0146 2.3000 0.0175 -0.1135 -0.1390 2.2500 -0.0127 -0.0506 -0.3358 1.9500 -0.0203 -0.1002 -0.5318 1.4000 1.2134 -0.1490 -0.7312 0.9000 -0.8393 0.9431 -0.4929 0.7000 0.0364 -0.0640 -0.1413 0.6000 -0.4480 0.0014 -0.1789 0.5000 0.5957 -0.5361 -0.3928 0.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*x20+.42770e-8*x19-.27708e-6*x18+.11098e-4*x17-.30784e-3*x16+.62785e-2*x15-.97558e-1*x14+1.1810*x13-11.297*x12+86.085*x11-524.68*x10+2558.0*x9-9942.3*x8+30586.*x7-73622.*x6+.13627e6*x5-.18907e6*x4+.18913e6*x3-.12803e6*x2+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; for n=1:N x(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-x2+1,x,y); t,u=Euler(ydot,0,0.5,0.025,20) 得到数值结果： t = Columns 1 through 17 0 0.0250 0.0500 0.0750 0.1000 0.1250 0.1500 0.1750 0.2000 0.2250 0.2500 0.2750 0.3000 0.3250 0.3500 0.3750 0.4000 Columns 18 through 21 0.4250 0.4500 0.4750 0.5000 u = Columns 1 through 17 0.5000 0.5375 0.5759 0.6153 0.6555 0.6966 0.7387 0.7816 0.8253 0.8700 0.9155 0.9618 1.0089 1.0569 1.1057 1.1553 1.2056 Columns 18 through 21 1.2568 1.3087 1.3613 1.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.4147 2. 改进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=1:N x(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-x2+1,x,y); t,u=Euler_ad(ydot,0,0.5,0.05,10) 得到数值结果： t = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 u = 0.5000 0.5768 0.6573 0.7414 0.8291 0.9202 1.0147 1.1126 1.2136 1.3178 1.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.4250 3.四阶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; for n=1:N x(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-x2+1,x,y); t,u=Runge_Kutta4(ydot,0,0.5,0.1,5) 得到数值结果： t = 0 0.1000 0.2000 0.3000 0.4000 0.5000 u = 0.5000 0.6574 0.8293 1.0151 1.2141 1.4256 结果比较： t Euler法 改进Euler法 四阶Runge_Kutta 精确解 0.1 0.6555 0.6573 0.6574 0.6574 0.2 0.8253 0.8291 0.8293 0.8293 0.3 1.0089 1.0147 1.0151 1.0151 0.4 1.2056 1.2136 1.2141 1.2141 0.5 1.4147 1.4250 1.4256 1.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(x3-x-1,3*x2-1) 当f(x)=x3-x-1; % x为初始点 % ep为精度，缺省值为1e-5 % it_max为最大迭代次数，缺省值100 % x_star为当迭代成功时输出方程的根，失败时输出最后的迭代值 % index为指标变量， index=1表示迭代成功 index=0表示失败 % it为迭代次数 if nargin4 it_max=100;end if nargin3 ep=1e-5;end index=0;k=1; while k=it_max x1=x; f=feval(fun,x); if abs(f(2)ep break; end x=x-f(1)/f(2); if abs(x-x1) fun=inline(atan(x),1/(1+x2); x_star,index,it=Newton(fun,1.45) 得到数值结果： x_star =1.6281e+004 index = 0 it =7 取初始值x0=0.5,输入 fun=inline(atan(x),1/(1+x2); x_star,index,it=Newton(fun,0.5) 得到数值结果： x_star =0 index = 1 it =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(x3-x-3,3*x2-1); x_star,index,it=Newton(fun,0.0) 得到数值结果： x_star = -0.0074 index = 0 it =101 取初始值x0=2.0,输入 fun=inline(x3-x-3,3*x2-1); x_star,index,it=Newton(fun,2.0) 得到数值结果： x_star = 1.6717 index = 1 it =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程序： 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 nargin4 it_max=100;end if nargin3 ep=1e-5;end n=length(A);k=0; x=zeros(n,1);y=zeros(n,1);index=1; while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end if abs(A(i,i)1e-10|k=it_max index=0;return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf) A,b=shuru;ep=1e-6; x,k,index=Jacobi(A,b,ep) 得到数值结果： x =(1.0000,1.0000,1.0000)T k = 19 index =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表示迭代失败。 if nargin5 it_max=100;end if nargin4 w=1;end if nargin3 ep=1e-5;end n=length(A);k=0; x=zeros(n,1);y=zeros(n,1);index=1; while 1 y=x; for i=1:n z=b(i); for j=1:n if j=i z=z-A(i,j)*x(j); end end if abs(A(i,i)1e-10|k=it_max index=0;return; end z=z/A(i,i);x(i)=(1-w)*x(i)+w*z; end if norm(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)T k = 11 index =1 依次取w=0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5 得到下表： 松弛因子w 迭代次数 0.8 19 0.9 16 1.0 13 1.1 11 1.2 14 1.3 17 1.4 23 1.5 27 结果分析： 从以上结果可以看出在求解相同问题时，sor法的迭代次数要比Jacobi迭代法少很多，计算量小许多。此外可以看出松弛因子w的选取对迭代次数的影响也十分大。在实际计算时，最优松弛因子很难事先确定，一般可用试算法取近似最优值。function X=jacobi(A,b,P,delta,max1) % A是n维非奇异阵 % B是n维向量 % P是初值 % delta是误差界 % X为所求的方程组AX=B的近似解 N=length(b); for k=1:max1 for j=1:N X(j)=(b(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j); end err=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 = 9 err = 2.5713e-005 X = 3.0000 2.0000 1.0000 nargin是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值，如下面的函数。例子，函数test1的功能是输出a和b的和。如果只输入一个变量，则认为另一个变量为0，如果两个变量都没有输入，则默认两者均为0。function y=test1(a,b)if nargin=0 a=0;b=0;elseif nargin=1 b=0;endy=a+b; s（2）迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。iJacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，即aii0(i=1,2,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角