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 - in

2、dex=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

3、 norm(y-x,inf)<ep break; end x=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.9998 11.9987 -3.0001k = 100index = 02.熟悉Gauss-Seidel迭代法，并编写Matlab程序function v,sN,vChain=gaussSeidel(A,b,x0,errorBound,maxSp)%Gauss-Seidel迭代法求解线性方程组%A-系数矩阵

4、b-右端向量 x0-初始迭代点 errorBound-近似精度 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,:)

5、=x0'k=k+1; error=norm(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.6169 11.1962 -4.2056sN = 11vChain = 6.0000 10.0000 -6.0000 -1.5000 2.0000 -3.5000 4.5000 10.3333 -5.5

6、000 -1.7500 3.6667 -3.4167 3.2500 10.6111 -5.0833 -1.9583 5.0556 -3.3472 2.2083 10.8426 -4.7361 -2.1319 6.2130 -3.2894 1.3403 11.0355 -4.4468 -2.2766 7.1775 -3.2411 0.6169 11.1962 -4.2056 0 0 0 0 0 0 0 0 0 0 0 0s数值实验 数值实验要求： 数值实验报告内容：要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。在本课程网站上提交数值实验报告的电子文档。 一、为了逼近飞

7、行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。 1 对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线； 2 对这些数据构造Lagrang插值多项式，并画出得到的Lagrang插值多项式曲线。 x 0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 f(x) 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 x 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3 f(x) 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.

8、5 0.4 0.25 1使用三次样条插值函数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 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

9、= 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.086

10、8 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

11、.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'); >

12、;> 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

13、,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.*

14、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的改进的

15、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=ze

16、ros(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

17、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

18、即采用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)

19、;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.1

20、000 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,

21、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

22、(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 精确解

23、 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迭代法求方程

24、 的根时，分别取初始值 ， ； 用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=in

25、line('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 nargin<4 it_max=100;end if nargin<3 ep=1e-5;end index=0;k=1; while k<=it_max x1=x; f=feval(fun,x); if abs(f(2)&l

26、t;ep break; end x=x-f(1)/f(2); if abs(x-x1)<ep index=1;break;end k=k+1; end x_star=x;it=k; 解：（1）由于f(x)=arctan(x) , f(x)= 1/1+x2 , 取初始值x0=1.45,输入 >> 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,输入 >

27、;> 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');

28、>> 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');

29、 得到图形： 结果分析： 从以上结果可以看出初值的选取与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

30、为最大迭代次数，缺省值100 % x为方程组的解 % k为迭代次数 % index为指标变量 index=1表示迭代收敛到指定要求 index=0表示迭代失败。 if nargin<4 it_max=100;end if nargin<3 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 i

31、ndex=0;return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf)<ep break; end x=y;k=k+1; end function A,b=shuru % 自己定义输入矩阵A和向量b的函数 n=15; for i=1:n if i=1|i=15 z(i)=3; else z(i)=2; end for j=1:n if j=i A(i,j)=4; elseif j=i+1|j=i-1 A(i,j)=-1; else A(i,j)=0; end end end b=z' 解：输入： >> A,b=shuru;e

32、p=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表示

33、迭代收敛到指定要求 index=0表示迭代失败。 if nargin<5 it_max=100;end if nargin<4 w=1;end if nargin<3 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

34、-w)*x(i)+w*z; end if norm(y-x,inf)<ep break; end k=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)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

35、 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);

36、 end err=abs(norm(X'-P); P=X' if (err<delta) break end end X=X'k,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是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值，如下面的函

37、数。 例子，函数test1的功能是输出a和b的和。如果只输入一个变量，则认为另一个变量为0，如果两个变量都没有输入，则默认两者均为0。 function y=test1(a,b) if nargin=0 a=0;b=0; elseif nargin=1 b=0; end y=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的下三角阵和上三角阵，于是Ax=b化为：x=D-1(L+U)x+D-1b与之对应