Jacobi迭代法 Gauss-Seidel迭代法
Matlab线性方程组的迭代解法（Jacobi迭代法Gauss-Seidel迭代法）实验报告2008年11月09日星期日12:49
1.熟悉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; 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
index=0; return;
end
y(i)=y(i)/A(i,i);
end
if norm(y-x,inf) 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.0001 k = 100 index = 2.熟悉Gauss-Seidel迭代法，并编写Matlab程序 function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp) %Gauss-Seidel迭代法求解线性方程组 %A-系数矩阵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 x0=inv(D)*(L+U)*x0+inv(D)*b; vChain(k,:)=x0'; k=k+1; error=norm(x0-fx0); fx0=x0; step=step+1; end v=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.2056 sN = 11 vChain = 6.0000 10.0000 -6.0000 -1.5000 2.0000 -3.5000 4.5000 10.3333 - 5.5000 -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 0 s 数值实验 数值实验要求： 数值实验报告内容：要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。在本课程网站上提交数值实验报告的电子文档。 一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。见下表。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.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 = 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'); 得到图形： https://www.360docs.net/doc/5018982873.html,grange插值算法： 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+.62785 e-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-9 942.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; 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-x^2+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-x^2+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-x^2+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('[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为迭代次数 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)) x=x-f(1)/f(2); if abs(x-x1) 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+x^2)]'); >> [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+x^2)]'); >> [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('[x^3-x-3,3*x^2-1]'); >> [x_star,index,it]=Newton(fun,0.0) 得到数值结果： x_star = -0.0074 index = 0 it =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.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程序： 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;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 index=0;return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf) 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;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 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-w)*x(i)+w*z; end if norm(y-x,inf) 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 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 nargin是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值，如下面的函数。 例子，函数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迭代法、超松弛迭代法和两步迭代法。 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迭代法的MA TLAB函数文件Jacobi.m如下： function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)>=eps x0=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迭代法的MA TLAB函数文件gauseidel.m如下： function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=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; %迭代次数 while norm(y-x0)>=eps x0=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]) 2013-2014（1）专业课程实践论文 题目：雅可比迭代法 一、算法理论 设有方程组),...,2,1(1 n i b x a i j n j ij ==∑= 记作,b Ax = （1） A 为非奇异阵且),,...,2,1(0n i a ij =≠将A 分裂为U L D A --=，其中 D =????????????????nn a a a 22 11，L =-??? ????? ???? ????-00001,21323121n n n n a a a a a a U =-?? ? ?? ? ? ? ????????-0000,122311312n n n n a a a a a a 将式（1）第)....2,1(n i i =个方程用ii a 去除再移项，得到等价方程组 (),,...,2,111n i x a b a x n i j j j ij i ii i =??? ? ? ?? -=∑≠= （2） 简记作 ,0f x B x += 其中 ().,111 0b D f U L D A D I B ---=+=-= 对方程组（2）应用迭代法，得到解式（1）的雅可比迭代公式 () () ()()()()()????????? ?? ? ??- ==∑≠=+,1,...,11002010n i j i k j ij i ii k i t n x a b a x x x x x ， 初始向量 （3） 其中()()()()()T k n k k k x x x x ,,...,21=为第k 次迭代向量。设()k x 已经算出，由式（3）可计算下一次迭代向量()(),,...,2,1,...;2,1,01n i k x k ==+ 显然迭代公式（3）的矩阵形式为 ()()()()???+=+,010f x B x x k k ，初始向量 其中0B 称为雅可比方法迭代矩阵。 第一题 算法解释 Jacobi迭代法 方程组Ax=b,其中A∈R nxn,b∈R n,且A为非奇异，则A可以写成A=D-L-U。 其中，D=diag[a11, a22,…, a nn],而-L，-U分别为A的上三角和下三角部分（不包括对角线元素） 则x=D-1(L+U)x+D-1b,由此可以构造迭代法：x(k+1)=Bx(k)+f 其中：B= D-1(L+U)x=I-D-1A，f=D-1b。 M文件 function[x,n]=jacobi(A,b,x0,eps,varargin) %采用Jacobi迭代法求线性方程组Ax=b的解 %线性方程组的系数矩阵：A %线性方程组中的常数向量：b %迭代初始向量：x0 %解的精度控制：eps %迭代步数控制：varargin %线性方程组的解：x %求出所需精度的解实际的迭代步数：n if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; n=1; %迭代次数 %迭代过程 while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end x n 例题 Jacobi迭代法求线性方程组实例。用Jacobi迭代法求解以下线性方程组10x1-x2=9 -x1+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 = 11 ans = 0.9958 0.9579 0.7916 实验五 矩阵的lu分解法，雅可比迭代法 班级： 学号： ： 实验五 矩阵的LU 分解法，雅可比迭代 一、目的与要求： ? 熟悉求解线性方程组的有关理论和方法； ? 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序； ? 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。 二、实验内容： ? 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解 各种方法的优缺点。 三、程序与实例 ? 列主元高斯消去法 算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+?表示 1) 消元过程 对k=1,2,…,n-1 ①选主元，找{}n ,,1k ,k i k +∈使得 k ,i k a = ik a n i k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。 ③如果k i k ≠，则交换第k 行与第k i 行对应元素位置， j i k j k a a ? j=k,┅,n+1 ④消元，对i=k+1, ┅,n 计算 kk ik ik a a l /= 对j=l+1, ┅,n+1计算 kj ik ij ij a l a a -= 2) 回代过程 ①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。 ②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算 ii n i j j ij n i i a x a a x /11,??? ? ? ?- =∑+=+ 程序与实例 程序设计如下： #include 软件迭代开发流程 前期项目引入，可行性分析略 项目调研：角色应包括项目经理、软件项目经理，应形成用户需求文档该文档需提交用户确认。产物为用户需求说明书文档 需求分析：角色应包括项目经理、软件项目经理、高级软件工程师，根据前期调研得到的用户需求说明书文档进行需求分析，应形成项目需求分析文档，该文档需提交项目组进行评审，主要是软件部，对需求能否实现进行评估。产物为项目需求分析说明书文档原型设计：角色应包括项目经理、UI设计、系统设计师，根据项目需求分析说明书进行原型设计，根据前期需求分析文档进行系统原型设计，主要包括利用界面原型制作工具设计图形类的功能模块，利用既有项目案例，制作实际项目案例参考，其中包括自己公司已有和市场上已存在的。连同项目需求分析说明书交由项目经理审核，最终由项目经理、软件项目经理同用户完成原型的审核，最终形成第一次迭代开发的项目需求文档说明书。 详细设计：角色应包括软件项目经理、项目组全体成员，应形成软件概要设计、软件详细设计文档，该文档需提交项目组，主要是项目部，对设计是否符合用户需求进行评估。经多次修改与确认后形成最终的项目详细设计说明书文档（包括概要设计）。产物为：项目概要设计说明书，项目详细设计说明书文档。 原型开发：角色应包括软件开发人员，按照详细设计说明书进行原型开发；快速实现详细设计说明中的各项功能，细节问题放到二次或三次迭代时加入。内部测试完毕后交由测试部进行测试。 测试评审：角色应为测试部、项目组成员，测试只进行功能实现测试，不进行其他细节和边界条件等测试。测试通过后，交由项目组进行评审，修改。最后由项目经理、软件项目经理与用户就原型进行沟通，检验功能设计是否符合用户要求，用户是否还有其他需求。最后形成二次迭代开发新需求文档，到此一次迭代结束。 流程图如下： 实验一非线性方程的数值解法（一）信息与计算科学金融崔振威201002034031 一、实验目的： 熟悉二分法和简单迭代法的算法实现。 二、实验内容： 教材P40 2.1.5 三、实验要求 1 根据实验内容编写二分法和简单迭代法的算法实现 2 简单比较分析两种算法的误差 3 试构造不同的迭代格式，分析比较其收敛性 (一)、二分法程序： function ef=bisect(fx,xa,xb,n,delta) % fx是由方程转化的关于x的函数，有fx=0。 % xa 解区间上限 % xb 解区间下限 % n 最多循环步数，防止死循环。 %delta 为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1:n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa) k=0; while abs(x-x0)>eps & k 作业六：分别编写用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组Ax=B的标准程序，并求下列方程组的解。 可取初始向量 X(0) =(0，0，0)’； 迭代终止条件||x(k+1)-x(k)||<=10e-6 （1） = （2） = Jacobi迭代法： 流程图 开 始 判断b中的最大值 有没有比误差大 给x赋初值 进行迭代 求出x，弱到100次还没到，警告不收 结束 程序 clear;clc; A=[8,-1,1;2,10,01;1,1,-5]; b=[1;4;3]; e=1e-6; x0=[0;0;0]'; n=length(A); x=zeros(n,1); k=0; r=max(abs(b)); while r>e for i=1:n d=A(i,i); if abs(d) 程序结果（1） （2） Gauss-Seidel迭代法： 程序 clear;clc; %A=[8,-1,1;2,10,01;1,1,-5]; %b=[1;4;3]; A=[5,2,1;-1,4,2;2,-3,10]; b=[-12;20;3]; m=size(A); if m(1)~=m(2) error('矩阵A不是方阵'); end n=length(b); %初始化 N=0;%迭代次数 L=zeros(n);%分解A=D+L+U,D是对角阵，L是下三角阵，U是上三角阵U=zeros(n); D=zeros(n); G=zeros(n);%G=-inv(D+L)*U d=zeros(n,1);%d=inv(D+L)*b x=zeros(n,1); for i=1:n%初始化L和U for j=1:n if i 数值计算方法实验报告（五） 班级：地信10801 序号：姓名： 一、实验题目：jacobi迭代法和Gauss-Seidel迭代法 二、实验学时: 2学时 三、实验目的和要求： 1．掌握迭代法的基础原理。 2．掌握jacobi迭代法和Gauss-Seidel迭代法的步骤。 3．能用程序语言对jacobi迭代法和Gauss-Seidel迭代法进行编程实现。 四、实验过程代码及结果 1、代码： #include for(int i=1;i<=N;i++) { float sum=0; for(int j=1;j<=N;j++) if(i!=j)sum+=xk[j]*a[i][j]; x[i]=(a[i][N+1]-sum)/a[i][i]; if(fabs(x[i]-xk[i]) printf("误差限是%f \n",tol); // 输出误差限 printf(" Jacobi迭代解序列X^(k) max|x^(k+1)-x^(k)| \n "); printf("x^%d = ",k=0); for(i=0;i 第六章 线性方程组的迭代法 一、教学目标及基本要求 通过对本节的学习，使学生掌握线性方程组的数值解法。 二、教学内容及学时分配 本节主要介绍线性方程组的数值解法，迭代公式的建立，迭代收敛性。 三、教学重点难点 1．教学重点：迭代公式的建立、迭代收敛性。 2. 教学难点：迭代收敛性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 6.2 解线性方程组的迭代法 重要性：解线性代数方程组的有效方法在计算数学和科学计算中具有特殊的地位和作用。如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题。在实际问题中产生的线性方程组的类型有很多，如按系数矩阵含零元素多少分类，有稠密和稀疏(零元素占80%以上)线性方程组之分；如按阶数的高低分类，有高阶(阶数在1000阶以上)中阶、(500~1000阶) 和低阶(500阶以下)线性方程组之分；如按系数矩阵的形状和性质分类，有对称正定、三对角、对角占优线性方程组之分。因为数值解法必须考虑方法的计算时间和空间效率以及算法的数值稳定性。因此，不同类型的线性方程组，其数值解法也不相同。但是，基本的方法可以归结为两大类，即直接法和迭代法。 分类：线性方程组的解法可分为直接法和迭代法两种方法。 (a) 直接法：对于给定的方程组，在没有舍入误差的假设下，能在预定的运算次数内求得精确解。最基本的直接法是Gauss 消去法，重要的直接法全都受到Gauss 消去法的启发。计算代价高。但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解，如何避免舍入误差的增长是设计直接法时必须考虑的问题。 (b) 迭代法：基于一定的递推格式，产生逼近方程组精确解的近似序列。收敛性是其为迭代法的前提，此外，存在收敛速度与误差估计问题。迭代法要求方程组系数矩阵具有某种特殊形式（如对角占优阵），是解高阶稀疏矩阵方程组的重要方法。 §6.1 迭代公式的建立 迭代法的基本思想是用逐次逼近的方法求线性方程组的解。 设有方程组b Ax = （1） 将其转化为等价的便于迭代的形式f Bx x += （2） （这种转化总能实现，如令b f A I B =-=,）并由此构造迭代公式 实验题目 1． Jacobi 迭代法 用Jacobi 迭代法求解线性方程组AX b =，保留四位有效数字(err ＝1e-4)，其中A=[8 -1 1;2 10 -1;1 1 -5]; b=[1 ;4; 3]。 2．Gauss-Seidel 迭代法 用Gauss-Seidel 迭代法求解线性方程组AX b =，保留四位有效数字(err ＝1e-4)，其中A=[8 -1 1;2 10 -1;1 1 -5]; b=[1 ;4; 3]。 3．分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组 123456410100014 10 105014001010041060101 4120 10 1 46x x x x x x ??--???? ??????---????????????--=??????--???????????? ----??????--?????? ?? 要求精度0.00001ε=，初始00=x ，最大迭代次数N=25，试比较这几种迭代法的迭代次数和收敛速度。 1.程序： #include MATLAB程序设计实践 1、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精 通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９ 年） “里查森迭代法线性方程组求解” 解： 算法说明： 里查森迭代法是最简单的迭代法，它的迭代公式为：x k+1=（I-A）*x k+b；在MATLAB 中编程实现的里查森迭代法函数为：richason。 功能：用里查森迭代法求线性方程组 调用格式：[x,n]=richason(A,b,x0,eps,M) 其中，A为线性方程组的系数矩阵； b为线性方程组的常数向量； x0为迭代初始向量； eps为解的精度控制（此参数可选）； M为迭代步数控制（此参数可选）； x为线性方程组的解； n为求出所需精度的解实际的迭代步数。 里查森迭代法的MA TLAB程序代码如下： function [x,n] = richason(A,b,x0,eps,M) %采用里查森迭代法求线性方程组Ax=b的解 %线性方程组的系数矩阵：A %线性方程组的常数向量：b %迭代初始向量：x0 %解的精度控制：eps %迭代步数控制：M %线性方程组的解：x %求出所需精度的解实际的迭代步数：n if(nargin==3) eps=1.0e-6; %eps表示迭代精度 M=200; %M表示迭代步数的限制值 elseif(nargin==4) M=200; end I=eye(size(A)); x1=x0; x=(I-A)*x0+b; n=1; %迭代过程 while(norm(x-x1)>eps) x1=x; x=(I-A)*x1+b; n=n+1; %n为最终求出解时的迭代步数 if(n>=M) 实验五 矩阵的lu分解法,雅可比迭代法 班级: 学号: 姓名: 实验五??矩阵的L Ｕ分解法,雅可比迭代 一、目的与要求: ? 熟悉求解线性方程组的有关理论和方法; ? 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序; ? 通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。 二、实验内容: ? 会编制列主元消去法、L Ｕ 分解法、雅可比及高斯—塞德尔迭代法德程序,进一步了解 各种方法的优缺点。 三、程序与实例 ? 列主元高斯消去法 算法：将方程用增广矩阵［Ａ∣b ］=（ij a )1n (n )+?表示 1) 消元过程 对k=1,2,…,n —１ ①选主元，找{}n ,,1k ,k i k +∈使得 k ,i k a =ik a n i k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③. ③如果k i k ≠,则交换第k 行与第k i 行对应元素位置, j i k j k a a ? j=k,┅，n+1 ④消元,对i=k+1, ┅,n 计算 kk ik ik a a l /= 对j=ｌ+1, ┅,n+1计算 kj ik ij ij a l a a -= 2) 回代过程 ①若0=nn a ，则矩阵Ａ奇异,程序结束;否则执行②。 ②nn n n n a a x /1,+=;对i=n —1, ┅,2,1，计算 ii n i j j ij n i i a x a a x /11,??? ? ??-=∑+=+ 程序与实例 程序设计如下: #include ＜iｏstreaｍ> ＃inｃlude <ｃmath> using ｎamespａce stｄ; voiｄ disｐ（ｄoｕｂle＊* p,ｉnt rｏｗ,inｔcｏl){ for（ｉnt i=0;i雅可比迭代法
Jacobi迭代法
雅可比迭代法与矩阵的特征值
软件迭代开发流程
二分法、简单迭代法的matlab代码实现教学文案
数值分析5-用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
jacobi迭代法和Gauss-Seidel迭代法
线性方程组的迭代法
实验4 Jacobi迭代法和GS迭代2
里查森迭代法线性方程组求解
雅可比迭代法与矩阵的特征值