fortran用高斯消去法矩阵求逆法、三角分解法、追赶法中的两种解线性方程组.doc

fortran:用高斯消去法，矩阵求逆法、三角分解法、追赶法中的两种解线性方程组 悬赏分：125 | 解决时间：2010-7-15 20:58 | 提问者：flygg888 一物理系统可用下列线性方程组来表示：从文件中读入m1、m2和的值，求a1、a2、N1 和N2的值。其中g取9.8，输入时以角度为单位。要求：（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三角分解法、追赶法等），定义求解线性方程组Ax=b的子程序，要求该子程序能求解任意线性方程组。（2）在主程序中分别调用上面定义的两个子程序，并对求解结果进行对比分析。（3）绘制以上两个方法所求得的方程解的数据分布图。请将答案发送至，答案合格者，将继续追加高分！问题补充： 请提供全面一点的答案（不要只求逆矩阵），谢谢了最佳答案 !高斯消去法 subroutine agaus(a,b,n,x,l,js) dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1 !逻辑变量 do k=1,n-1 d=0.0 do i=k,n do j=k,n if (abs(a(i,j)d) then d=abs(a(i,j) js(k)=j is=i end if end do end do !把行绝对值最大的元素换到主元位置 if (d+1.0=1.0) then l=0 else !最大元素为0无解 if(js(k)/=k) then do i=1,n t=a(i,k) a(i,k)=a(i,js(k) a(i,js(k)=t end do !最大元素不在K行，K行 end if if(is/=k) then do j=k,n t=a(k,j) a(k,j)=a(is,j) a(is,j)=t !交换到K列 end do t=b(k) b(k)=b(is) b(is)=t end if !最大元素在主对角线上 end if !消去 if (l=0) then write(*,100) return end if do j=k+1,n a(k,j)=a(k,j)/a(k,k) end do b(k)=b(k)/a(k,k) !求三角矩阵 do i=k+1,n do j=k+1,n a(i,j)=a(i,j)-a(i,k)*a(k,j) end do b(i)=b(i)-a(i,k)*b(k) end do end do if (abs(a(n,n)+1.0=1.0) then l=0 write(*,100) return end if x(n)=b(n)/a(n,n) do i=n-1,1,-1 t=0.0 do j=i+1,n t=t+a(i,j)*x(j) end do x(i)=b(i)-t end do100 format(1x,fail) js(n)=n do k=n,1,-1 if (js(k)/=k) then t=x(k) x(k)=x(js(k) x(js(k)=t end if end do return end program maindimension a(4,4),b(4),x(4),js(4)double precision a,b,xreal m1,m2,jopen(1,file=laiyi.txt)read(1,*)m1,m2,jclose(1) n=4print*,m1,m2,ja(1,1)=m1*cos(3.14159*j/180)a(1,2)=-m1a(1,3)=-sin(3.14159*j/180)a(1,4)=0a(2,1)=m1*sin(3.14159*j/180)a(2,2)=0a(2,3)=cos(3.14159*j/180)a(2,4)=0a(3,1)=0a(3,2)=m2a(3,3)=-sin(3.14159*j/180)a(3,4)=0a(4,1)=0a(4,2)=0a(4,3)=-cos(3.14159*j/180)a(4,4)=1b(1)=0 b(2)=m1*9.8b(3)=0b(4)=m2*9.8 call agaus(a,b,n,x,l,js) if (l/=0) then write(*,*)a1=,x(1),a2=,x(2) ,n1=,x(3),n2=,x(4)end ifend!逆矩阵求解SUBROUTINE qiuni(A,N,L,IS,JS)DIMENSION A(N,N),IS(N),JS(N)DOUBLE PRECISION A,T,DL=1DO K=1,N D=0.0 DO I=K,N DO J=K,N IF(ABS(A(I,J).GT.D) THEN !把最大的元素给D D=ABS(A(I,J) IS(K)=I JS(K)=J END IF END DO END DO IF (D+1.0.EQ.1.0)THEN L=0 WRITE(*,200) RETURN END IF 200 FORMAT(1X,ERR*NOT INV) DO J=1,N T=A(K,J) A(K,J)=A(IS(K),J) A(IS(K),J)=T END DO DO I=1,N T=A(I,K) A(I,K)=A(I,JS(K) A(I,JS(K)=T END DO A(K,K)=1/A(K,K) DO J=1,N IF(J.NE.K)THEN A(K,J)=A(K,J)*A(K,K) END IF END DO DO I=1,N IF(I.NE.K)THEN DO J=1,N IF(J.NE.K)THEN A(I,J)=A(I,J)-A(I,K)*A(K,J) END IF END DO END IF END DO DO I=1,N IF(I.NE.K)THEN A(I,K)=-A(I,K)*A(K,K) END IF END DOEND DODO K=N,1,-1 DO J=1,N T=A(K,J) A(K,J)=A(JS(K),J) A(JS(K),J)=T END DO DO I=1,N T=A(I,K) A(I,K)=A(I,IS(K) A(I,IS(K)=T END DO END DO RETURN ENDSUBROUTINE BRMUL(A,B,N,C)DIMENSION A(N,N),B(N),C(N)DOUBLE PRECISION A,B,CDO I=1,NDO J=1,NC(I)=0.0DO L=1,NC(I)=C(I)+A(I,L)*B(L)END DO END DOEND DORETURNENDprogram mainDIMENSION A(4,4),B(4,1),C(4,1),IS(4),JS(4)DOUBLE PRECISION A,B,CREAL M1,M2,JDOPEN(1,FILE=LAIYI.TXT)READ(1,*) M1,M2,JDPRINT*,M1,M2,JDCLOSE(1)A(1,1)=M1*COS(3.14*JD/180)A(1,2)=-M1A(1,3)=-SIN(3.14*JD/180)A(1,4)=0A(2,1)=M1*SIN(3.14*JD/180)A(2,2)=0A(2,3)=COS(3.14*JD/180)A(2,4)=0A(3,1)=0A(3,2)=M2A(3,3)=-SIN(3.14*JD/180)A(3,4)=0A(4,1)=0A(4,2)=0A(4,3)=-COS(3.14*JD/180)A(4,4)=1B(1,1)=0B(2,1)=M1*9.8B(3,1)=0B(4,1)=M2*9.8CALL QIUNI(A,4,L,IS,JS)CALL BRMUL(A,B,4,C)WRITE(*,*) (C(I,1),I=1,4)END画图 USE MSFLIB INTEGER status TYPE(xycoord) xy status=SETCOLORRGB(#FFFFFF) status1=SETCOLORRGB(#0000FF) OPEN(1,FILE=G.TXT) READ(1,*) G1,G2,G3,G4 OPEN(2,FILE=N.TXT) READ(2,*) N1,N2,N3,N4 CALL MOVETO(INT(20),INT(20),xy) status=LINETO(INT(40),INT(G1) status=LINETO(INT(80),INT(G2) status=LINETO(INT(120),INT(G3) status=LINETO(INT(160),INT(G4) CALL SE