Jacobi迭代法.docx

1、Jacobi迭代法第一题算法解释Jacobi迭代法方程组Ax=b，其中ARnxn，bRn，且A为非奇异，则A可以写成A=D-L-U。其中，D=diaga11, a22, , ann，而L，U分别为A的上三角和下三角部分(不包括对角线元素）则x=D1(L+U)x+D-1b，由此可以构造迭代法:x(k+1)=Bx（k)+f其中：B= D-1(L+U）x=I-D-1A，f=D-1b.M文件functionx,n=jacobi(A，b，x0,eps,varargin)%采用Jacobi迭代法求线性方程组Ax=b的解线性方程组的系数矩阵：A%线性方程组中的常数向量：b迭代初始向量：x0%解的精度控制：e

2、ps%迭代步数控制:varargin%线性方程组的解：x%求出所需精度的解实际的迭代步数:nif nargin=3 eps=1.0e6; M=200;elseif nargin3 error returnelseif nargin=5 M=varargin1;endD=diag(diag(A）); %求A的对角矩阵L=-tril(A，1); 求A的下三角阵U=-triu（A，1）； %求A的上三角阵B=D(L+U）;f=Db;x=B*x0+f；n=1； 迭代次数迭代过程while norm（xx0）=eps x0=x; x=Bx0+f; n=n+1; if(n=M) disp（Warning:

3、迭代次数太多，可能不收敛！）； return； endendxn例题Jacobi迭代法求线性方程组实例。用Jacobi迭代法求解以下线性方程组10x1x2=9x1+10x2+2x3=7-x2+10x3=6 在matlab命令窗口输入如下程序： a=10 -1 0；1 10 -2;0 2 10； b=9;7；6； jacobi（a，b，0;0；0）y = 0.9958 0.9579 0。7916输出的迭代次数为：n = 11ans = 0.9958 0.9579 0.7916n=length(x） 第二题方法一：运用戴维宁定理求出外电路的等效电压Uoc=(Z2/(Z1+Z2）-Z4/(Z3+Z4

4、）).Us，将电压源短路，z1,z2并联，z3，z4并联，Req=Z3。Z4./(Z3+Z4)+Z1。*Z2./（Z1+Z2）。所求电流源两端电压为Ut=Is。Req+Uocw=eps，1，2；Us=10,10,0;Is=5，5,0；Z1=1。/(0.5*wj）;Z4=1*w*j;Z2=2,2，2；Z3=2,2,2;Uoc=（Z2/（Z1+Z2)-Z4/(Z3+Z4)）.*Us；Req=Z3.Z4./（Z3+Z4)+Z1.*Z2./(Z1+Z2）;Ut=Is。Req+Uocdisp( w Um phi )disp（w，abs(u),angle(u)180/pi)等效电路图方法二:运用电路叠加定

5、理，分别求出电流源与电压源单独作用在所求端点两端的电压，当只有电压源作用时，电流源断路，得到Uoc=(z2./(z1+z2）z4。/(z3+z4）。us;当只有电流源作用时，电压源短路,z1,z2并联,z3，z4并联,外加电阻（z3。*z4./(z3+z4)+z1.z2。/(z1+z2），电流源作用得到Uc=Is.*。所求端点电压为两者叠加值Ut w=eps，1，2；输入频率值 us=10，10，0;%输入不同频率下的电压分量 Is=5，0,5； z1=1./（0。5w*j); z4=1w*j; z2=2,2,2; z3=2，2,2; Uoc=（z2。/(z1+z2）-z4。/(z3+z4)）

6、。us；%电压源单独作用的端点电压 Uc=Is.(z3.*z4。/(z3+z4）+z1.z2./(z1+z2）)； 电流源单独作用的端点电压 Ut=Uoc+Uc； disp(w，abs(u)，angle（u)180/pi） 0。0000 10.0000 0 1.0000 3。1623 18。43492.0000 7.0711 -8.1301开始第三题:已知f(x）=3sin(pi*x/6).2,x=1。5，3.5 k Xk f(Xk） 0 1 2 3 4 0。0 1.0 2.0 3。0 4.0 0.000 0。75 2.25 3.0 2.25求：(1)计算函数差商表。 （2）写出牛顿多项式Nk

7、（x),k=1，2，3,4. （3)在给定值x处求牛顿多项式Nk(x）,k=1,2,3,4的值。 （4)比较（3）中的结果与实际函数值。M文件functionP,A=chashang(X，Y） n=length（X）；A=zeros(n,n）； A（:，1）=Y; for j=2:n for i=j:n A(i，j）=（A(i,j-1)A（i1，j-1)）/（X（i）X（ij+1）); end endP=A(1:n,1：n)在命令窗口中输入 X=0 1 2 3 4; Y=0 0.75 2。25 3.0 2.25; chashang（X，Y）P = 0 0 0 0 0 0。7500 0。7500

8、 0 0 0 2。2500 1。5000 0。3750 0 0 3.0000 0。7500 -0.3750 -0.2500 0 2。2500 -0.7500 0.7500 0。1250 0。0313ans = 0 0 0 0 0 0.7500 0。7500 0 0 0 2。2500 1.5000 0.3750 0 0 3.0000 0.7500 0。3750 -0。2500 0 2.2500 0.7500 -0.7500 0。1250 0.0313差商程序流程图 M文件function C,A=newtonpoly(X，Y) n=length（X）；A=zeros（n,n)； A（：，1）=Y

9、； for j=2:n for i=j：n A(i,j）=（A（i，j-1）-A（i1,j1)）/(X（i)-X(ij+1）； end endC=A（n,n);for k=(n1):1：1 C=conv(C,poly（X(k）); d=length(C)； C(d)=C(d）+A（k，k)；end在命令窗口中输入 X=0 1 ； Y=0 0。75； n1=newtonpoly(X，Y)n1 = 0。7500 0 N1=poly2sym(n1)N1 =3/4*x X=0 1 2； Y=0 0.75 2。25； n2=newtonpoly(X，Y)n2 = 0。3750 0。3750 0 N2=p

10、oly2sym（n2）N2 =3/8x2+3/8*x X=0 1 2 3； Y=0 0.75 2.25 3.0； n3=newtonpoly（X，Y)n3 = 0。2500 1.1250 0。1250 0 N3=poly2sym(n3)N3 =1/4*x3+9/8x2-1/8*x X=0 1 2 3 4; Y=0 0。75 2。25 3。0 2.25; n4=newtonpoly(X，Y）n4 = 0.0313 -0。4375 1.4688 -0。3125 0 N4=poly2sym（n4）N4 =1/32*x4-7/16x3+47/32x25/16x f1=polyval(n1,1.5)f1

11、 = 1.1250 f2=polyval(n2，1。5)f2 = 1.4063 f3=polyval（n3,1.5）f3 = 1。5000 f4=polyval（n4,1。5)f4 = 1.5176 p1=polyval（n1,3.5)p1 = 2。6250 p2=polyval（n2,3.5）p2 = 5.9063 p3=polyval（n3，3。5)p3 = 2.6250 p4=polyval(n4，3.5）p4 = 2.8301 x=1。5； F=3*（sin(pix/6）)2F = 1.5000 x=3.5； P=3*（sin(pix/6）)2P =2.7990画图比较输入:x=0：0。05:4；y=3.*sin（pi.*x。/6)。sin（pi.*x./6); plot(x,y，r) hold on x=1.5； f1=1。125； f2=1.4063 f3=1。500; f4=1.5176 plot(x，f1，o，x,f2，o,x，f3，o,x,f4，o） x=3。5； p1=2。6250； p2=5。9063； p3=2。6250; p4=2.8301； plot（x，p1，o,x,p2,o，x，p3,o,x，p4,o）得到如下图形