矩阵运算2感应测井的半解析半数值计算方法.ppt

* * §4.3线性代数方程组的求解 所以求解过程分为两步: (1)消元 (2)回代 为了提高计算的精度，要采用选主元的技巧，有列主元法、行主元法和全主元法。 列主元高斯消去法程序见下页 GS.FOR SUBROUTINE GS(N,M,A,B) IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(N,N),B(N,M),P,T DO 70 K=1,N P0=0.D0 DO 10 I=K,N P1=ABS(A(I,K)) IF(P1.LE.P0)GOTO 10 P0=P1 I0=I 10 CONTINUE IF(P0.LE.1.D-16)PAUSE 9999 IF(I0.EQ.K)GOTO 40 DO 20 J=K,N T=A(K,J) A(K,J)=A(I0,J) 20 A(I0,J)=T DO 30 IM=1,M T=B(K,IM) B(K,IM)=B(I0,IM) 30 B(I0,IM)=T 40 P=1/A(K,K) K1=K+1 DO 50 J=K1,N A(K,J)=A(K,J)*P DO 50 I=K1,N 50 A(I,J)=A(I,J)-A(I,K)*A(K,J) DO 60 I=1,M 60 B(K,I)=B(K,I)*P DO 70 J=K1,N DO 70 I=1,M 70 B(J,I)=B(J,I)-B(K,I)*A(J,K) DO 80 I=N-1,1,-1 DO 80 J=I+1,N DO 80 K=1,M 80 B(I,K)=B(I,K)-A(I,J)*B(J,K) END ?若方程组的左端矩阵具有正定对称性，采用克劳特法效率最高，因为它只有约N 3/6的乘法次数。若矩阵是带形的，可进一步减少计算量。求解线性代数方程组的方法共有十几种。可针对不同的情况，选用不同的求解方法。 ?若方程组具有对称性，还可减少一半的计算量。其它求解方法还有 Gauss-Jordan 消去法,乔累斯基法、克劳特法等许多方法。其实，那些方法都是高斯消去法的变种。 这里我们还要介绍一个基于矩阵三角分解的求解线性方程组的方法. 如果 A 是一个上三角矩阵或下三角矩阵,我们可以用简单的回代公式求出方程的解. 对一般矩阵 A, 如果它能分解成一个上三角矩阵 U 和一个下三角矩阵 L 的乘积, 方程也可容易求出. 该方法称为 LU 分解法 可以证明一般矩阵 A 可分解成 其中 U 是单位上三角矩阵(对角元素为1). L,U的元素可用下列递推公式求得 线性方程可写成 回代公式 最后要注意病态线性方程组 良态方程 其精确解是 若方程系数有一点小扰动 其精确解是 病态方程 其精确解是 若方程系数有一点小扰动 其精确解是 如何判断方程 Ax=B 的性态? 可以根据系数矩阵的条件数来判断 条件数越大, 方程的病态越严重 对上面的良态方程的例子 对上面的病态方程的例子 §4.4感应测井的半解析半数值计算方法 图4.1地层结构和测量线圈剖面示意图 R井眼 sm T sx4 sx3 sx2 sx1 st1 st2 st3 st4 一、物理模型。 我们假设： a. 分界面是完全水平的； b. 分界面两侧的电导率是突变的； c.井壁和侵入带界面完全是同轴圆筒形，且与水平界面完全垂直； d.在每一个子区域内，电导率是完全均匀的； e.电极棒的中轴与井眼中轴重合，并忽略电极棒的存在； f.线圈是圆形的，并与中轴完全垂直，线圈没有宽度和厚度。因此，整个问题具有旋转对称性。 这就是我们对该问题所作的物理模型 图4.1地层结构和测量线圈剖面示意图 R井眼 sm T sx4 sx3 sx2 sx1 st1 st2 st3 st4 二、数学模型。 上述方程和边界条件就是该问题的数学模型。 三、计算方法。 该问题一般没有解析解。只有两个特殊情况才有解析解：无井无侵水平多分层情况和有井有侵无限厚单地层同轴圆柱面情况。这两种情况下的解析解是用无穷积分来表达的。只有单一均匀介质情况下的解才是一个有限的表达式。 一种数值方法是有限元法，80年代用大型计算机算需时25小时。 另一种数值方法是半解析半数值方法。实际计算的结果表明，计算穿过10米厚地层的问题，在386-33微机上仅需时半分钟。 a.分离变量得部分解析解。 b.齐次方程的通解。