Matlab第一次上机作业
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
输入: >>tic, n=9;[u,k]=xsj(n), toc,surf(u)
计算Байду номын сангаас果如下:
u=
0 0 0 0 0 0 0 0 0 0 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0577 0.0986 0.1258 0.1412 0.1462 0.1412 0.1258 0.0986 0.0577 0 0 0.0560 0.0955 0.1216 0.1364 0.1412 0.1364 0.1216 0.0955 0.0560 0 0 0.0508 0.0859 0.1088 0.1216 0.1258 0.1216 0.1088 0.0859 0.0508 0 0 0.0413 0.0686 0.0859 0.0955 0.0986 0.0955 0.0859 0.0686 0.0413 0 0 0.0256 0.0413 0.0508 0.0560 0.0577 0.0560 0.0508 0.0413 0.0256 0 0 0 0 0 0 0 0 0 0 0 0 0
程序一: function [u,k]=xsbj(n) % xsbj:用块 Jacobi 迭代法求解线性方程组 A*u=f % u:方程组的解; k 迭代次数; n:非边界点数 % a:方程组系数矩阵 Aii 的下对角线元素;b:方程系数矩阵 Aii 的主对角线元素 % c:方程组系数矩阵 Aii 的上对角线元素;d:追赶法所求方程的右端向量 % e:允许误差界;er:迭代误差 f=2*1/(n+1)^(2)*ones(n+2,n+2); a=-1*ones(1,n);b=4*ones(1,n);c=-1*ones(1,n);u=zeros(n+2,n+2);e=0.000000001; for k=1:2000 er=0; ub=u; for j=2:n+1 d(1:n)=f(2:n+1,j)+ub(2:n+1,j-1)+ub(2:n+1,j+1); x=zg(a,b,c,d); u(2:n+1,j)=x'; er=er+norm(ub(:,j)-u(:,j),1); end if er/n^2<e,break; end end 程序二: function x=zg(a,b,c,d) % zg:用追赶法求解线性方程组 n=length(b); % LU 分解 u(1)=b(1); for k=2:n if u(k-1)==0,D=0,return; end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); end % 追赶法求解之追过程,求解 Ly=d y(1)=d(1); for k=2:n y(k)=d(k)-l(k)*y(k-1); end % 追赶法求解之追过程,求解 Ux=y if u(n)==0,D=0,return;end x(n)=y(n)/u(n); for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1))/u(k); end 输入:
将差分方程组写成矩阵形式为������������ = ������,其中������ ∈ ������
������× ������ × ������× ������
, ������, ������ ∈ ������ ������× ������ 。
差分方程等价于 Av=b。这样,原问题就可用 Jacobi 迭代法或块 Jacobi 迭代法求解方程 组。块 Jacobi 迭代法公式可化简为 A������������ ������������������ +1 = ������������
用 MATLAB 语言编写求解线性方程组 Au=f 的算法程序, 采用下列方法计算结 果和算法性能,对计算结果给出结论。 (1)用 Jacobi 迭代法求解线性方程组 Au=f; (2)用块 Jacobi 迭代法求解线性方程组 Au=f; (3)用(预条件)共轭斜量法求解线性方程组 Au=f。
解: 1、基本原理说明:
(������ )
+ ������������ +1
(������ )
+ ������������
接下来就可以编写适当的数值算法求解线性方程组。
2、计算方法:
(1)Jacobi 迭代法:算法程序如下: function [u,k]=xsj(n) % xsj:用 Jacobi 迭代法求解线性方程组 A*u=f % u:方程组的解; k 迭代次数; n:非边界点数 % e:允许误差界;er:迭代误差 f(2:n+1,2:n+1)=(n+1)^(-2)*2; u=zeros(n+2,n+2); e=0.000000001; for k=1:1000 er=0; ub=u; for j=2:n+1 for i=2:n+1 u(i,j)=(ub(i-1,j)+ub(i+1,j)+ub(i,j-1)+ub(i,j+1)+f(i,j))/4; er=er+abs(ub(i,j)-u(i,j)); end end if er/n^2<e,break; end end
首先剖分求解区域。用平行于坐标轴的直线 ������ = ������������ = ������ℎ, ������ = ������������ = ������ℎ, ℎ = 1 ������ + 1
������, ������ = 0,1,· · · , ������, ������ + 1 将求解区域分为网格。然后用数值微分公式对微分算子进行离散。 ������ 2 ������ ������ ������������−1 , ������������ − 2������ ������������ , ������������ + ������ ������������ +1 , ������������ ������ , ������ = + ������ ℎ2 ������������ 2 ������ ������ ℎ2 ������ 2 ������ ������ ������������ , ������������−1 − 2������ ������������ , ������������ + ������ ������������ , ������������ +1 ������������ , ������������ = + ������ ℎ2 2 ������������ ℎ2 即有 ������������−1,������ − 2������������ ,������ + ������������ +1,������ ������ 2 ������ ������������ , ������������ ≈ 2 ������������ ℎ2 ������������ ,������ −1 − 2������������ ,������ + ������������ ,������ +1 ������ 2 ������ ������ , ������ ≈ ������������ 2 ������ ������ ℎ2 其中,uij 表示 u(xi, yi)的近似值。得到在每个点(xi, yi)上的有限差分方程为: 4������������ ,������ − ������������−1,������ − ������������ +1,������ − ������������ ,������ −1 − ������������ ,������ +1 = ℎ2 ������������������ 1 ≪ ������ , ������ ≪ ������ (3-74) 又称五点差分格式。其中 fij=f(xi, yi)。在边界上,有 ������0,������ = ������������ +1,������ = ������������ ,0 = ������������ ,������+1 = 0 1 ≪ ������ , ������ ≪ ������ (3-75)
k =304 Elapsed time is 0.005737 seconds.
0.2
0.15
0.1
0.05
0 15 10 10 5 0 0 5 15
Figure 1 Jacobi 迭代法数值解图示
Jacobi 迭代法最终迭代结果与教材中 Gauss-Seidel 迭代法结果一致,但是 Jacobi 法的迭 代次数为 304,比 Gauss-Seidel 迭代法要大近一倍,而且 Gauss-Seidel 法在本人计算机上运 行时间为 0.003540(多次平均,下同) ,Jacobi 法为 0.005737。可见 Gauss-Seidel 迭代法的动 态存储方式不仅节省了一半的存储空间, 而且加快了收敛速度, 比 Jacobi 法更早的得出结果。 (2)块 Jacobi 迭代法:算法程序如下:
对非边界点编号,顺序为对直线 y=yi 从上往下,对固定的一条直线的网点从左往右依 次编号,即有 ������1 , ������1 , ������2 , ������1 ,· · · , ������������ , ������1 , ������1 , ������2 , ������2 , ������2 , ������2 , ������2 ,· · · , ������������ , ������2 ,· · · , ������1 , ������������ , ������2 , ������������ ,· · · , ������������ , ������������ 相应的解向量和右端向量分别为 ������ = (������1,1 , ������2,1 ,· · · , ������������ ,1 , ������1,2 , ������2,2 ,· · · , ������������ ,2 ,· · · , ������1,������ , ������2,������ ,· · · , ������������ ,������ )������ ������ = ℎ2 (������1,1 , ������2,1 ,· · · , ������������ ,1 , ������1,2 , ������2,2 ,· · · , ������������ ,2 ,· · · , ������1,������ , ������2,������ ,· · · , ������������ ,������ )������ 由 3-74 和 3-75 得差分方程组 4������������ ,������ − ������������−1,������ − ������������ +1,������ − ������������ ,������ −1 − ������������ ,������ +1 = ℎ2 ������������������ ������0,������ = ������������ +1,������ = ������������ ,0 = ������������ ,������+1 = 0 1 ≪ ������ , ������ ≪ ������
第一次上机作业
数值实验题三 2.用有限差分法(五点分格式)求解正方形域上的 Poission 方程边值问题
− ������������ ������ ������������ ������ + = ������ ������, ������ = ������ ������������������ ������������������ ������ ������, ������ = ������ ������, ������ = ������ ������, ������ = ������ ������, ������ = ������