在计算机上采用动态计算形式
??(k?1)(k?1)(k)(k)(k?1) du?(bi,j?ui?1,j?ui,j?1?ui,j?1?ui?1,j?4ui,j)?4?
?ui,j?ui,j?du(i,j?2，3，???，n)?
如果|du|>error则error=|du|,如果error 第三步：比较不同的w的迭代次数，用kk存放最小迭代次数，用ww和uu存放相应 的w及u。 3． 程序 ① [u,k]=SOR(u,b,w) %%%%%%%（被下面程序调用） %输入场初值u0、场源b及松弛因子w，通过五点差分格式进行迭代运算， %如果第k+1次的迭代结果与第k次的差小于精度，则可以近似认为第k+1次的迭代 %结果是精确解，然后返回迭代次数k和迭代解 function [u,k]=SOR(u,b,w) %输出迭代结果u，及迭代次数k m=length(u); %m为u的维数 h=1/(m-1); %h为步长 N=10000;e=0.0000001; %e为精度 for k=1:N %k为记录迭代次数 error=0; for j=2:m-1 for jj=2:m-1 sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1); du=w*(h^2*b(jj,j)-sum)/4; %计算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if error end if error<=e,break;end %判断是否达到精度 end ② [kk,ww,uu]=SOR_5dianchafen(n) %用超松弛迭代法求解正方形域上的Poisson方程边值问题 %用5点差分格式求取泊松问题 %输入n,对x、y轴进行n等分；先确定场u的边界及场源b,在调用[u,k]=SOR(u,b,w)； %用不同w计算的迭代次数不同，用kk存放最小的迭代次数， %用ww和uu分别存放最佳松弛因子w和精确解 function [kk,ww,uu]=SOR_5dianchafen(n) w=[1.1:0.1:1.8];m=length(w); %w为松弛因子 kk=1000; ww=0; %kk是最少迭代次数，ww是最松弛因子 h=1/n; %h步长 b=zeros(n+1,n+1); %计算场源b tic; for i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h^2; end end uu=zeros(n+1,n+1); u=zeros(n+1,n+1); %对u赋初值 u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1; mu=u; %初值mu以便不同的w计算 for i=1:m %用不同的w计算迭代 [u,k]=SOR(mu,b,w(i)); %调用[u,k]=SOR(u,b,w)，返回迭代次数及精确解 if kk>k, kk=k;ww=w(i);uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uu end t=toc %统计程序运算时间 4.计算结果： >> format short >> n=10; >> [kk,ww,uu]=SOR_5dianchafen(n) t = 0.0310 kk = 48 ww = 1.6000 uu = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 >> contourf (uu, 'DisplayName', 'uu'); figure(gcf) 11109876543211234567891011 图一 超松弛 二 用块Jacobi迭代法求解线性方程组Au=f。 1. 基本原理： 对A做自然分解A=D-L-U=D-(L+U) 其中D是有A的对角线元素组成的矩阵，L是由A的对角线以下元素组成的矩阵，U是由A得对角线以上元素组成的矩阵。 于是将M=D,N=L+U,代入得到 Dx=(L+U)x+b 任取x的初值进行迭代 2. 算法： （1）Gauss-Sediel迭代法原理 i?1 (k?1)(k?1) xi?(bi?n?aj?1xjij???uaijxj)/aii,(k)j?i?1i?1,2,?,n五点差分格式： (k?1)(k?1)(k?1)(k)(k)i,ji?1,ji,j?1i,j?1i?1,j i,j因为A可以写成块状，即： u?(b?u?u?u)/4?A212??IA???????IA33????I????I??AN?1,N?1??T 如果把每一条线上的节点看作一个组v3?(u2,3,u3,3,...,uN,3)，可以把Au=f表示成块状求解： Ajjvj?vj?1?vj?1?bj, （2）计算步骤： 2?j?N 第一步：给场值u和场源b赋初值，及定义 第二步：用公式 Ajjvj?vj?1?vj?1?bjAii ，进行迭代计算 第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。 3. 程序 [k,u]=kuai_GaussSeidel(n) %用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题 %输入n,对x、y轴进行n等分；先确定场u的边界及场源b； %用k和u分别存放迭代次数和精确解 function [k,u]=kuai_GaussSeidel(n) %对x、y轴进行n等分 h=1/n; %步长 u=zeros(n+1,n+1); %对u赋初值 u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1; b=zeros(n+1,n+1); %计算场源b for i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h^2; end end b=h^2*b; for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1); b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1); end A=zeros(n-1,n-1); %定义矩阵的子块A for i=1:n-1 if i>1, A(i,i-1)=-1; end if i e=0.000000001; %e是精度 tic; for k=1:1000 %k是迭代次数 error=0; %误差 for j=2:n %进行迭代循环 ub=u; %ub是第k-1次迭代结果，用于和第k次迭代结果比较 if j==2, u(2:n,2)=pinv(A)*(u(2:n,3)+b(2:n,2)); end if j==n, u(2:n,n)=pinv(A)*(u(2:n,n-1)+b(2:n,n)); end if j>2&j error=max(max(abs(u-ub))); %error是前后两次迭代结果的对应元素的最大误差 if error<=e, break; end %判断误差是否达到精度 end t=toc %统计程序运算时间 4. 计算结果： >> format short >> n=10; >> [k,u]=kuai_GaussSeidel(n) t = 1.3280 k = 93 u = Columns 1 through 8