结构力学课程设计专业：班级：姓名：学号：指导老师：日期：2015年7月5日目录前言 (1)问题一： (3)问题描述： (3)程序说明： (3)全选主元高斯约当消去法： (3)全选主元高斯约当消去法的程序及注解如下： (4)运行结果： (6)问题二： (6)问题描述： (6)方法一：追赶法 (7)程序说明： (7)追赶法带型的计算程序及注解： (7)运行结果： (9)总结与思考： (9)方法二：列选主元高斯消去法算带型问题 (10)程序说明： (10)列选主元高斯消去法算带型计算程序及注解： (10)运行结果： (12)反思与对比（收获）： (12)问题三： (13)问题描述： (13)程序框图： (14)程序特点： (14)1.主要变量： (15)2.子例行子程序哑元信息： (15)3.文件管理： (16)4.数据文件格式： (16)源程序： (17)输入数据如下（input.txt）： (23)输出数据如下（output.txt）： (23)程序运行后输出数据结果如下（需要手动打开output.txt文件）： (24)总结与收获： (25)参考文献: (26)前言：经过这学期的学习与积累，对结构力学这门课程有所收获，结构力学这门课程对我们学习飞行器设计与专业的学生来说，那就是手足的关系，因为我感觉任何航空、航天器都离不开结构的设计，只要有结构就牵涉到结构力学的分析与计算，因为航空器在空中飞行要遇到很多“挫折”，结构力学就是来分析这些个“挫折”下，看航空器能不能经受得了。

结构力学课程从内容上讲，主要涉及机构的几何组成分析，求解静定、超静定结构内力的虚功原理。

具体分析问题的方法包括力法、位移法等。

但对于复杂结构来讲，简单的手算的方法过于繁琐。

因此，由于课程设计偏重于利用Fortran 语言编写有限元子程序来完成复杂结构的内力计算，我就恶补了好几天的与Fortran有关的知识，下面就现学现卖的计算了王老师给的三个问题，肯定有不妥之处，希望读者纠错。

问题一：一、 利用全选主元的高斯约当（ Gauss-Joadan ）消去法求解如下方程组，并给出详细 的程序注解和说明：⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧=⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧∙⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡1536353424543214109753910862781071567554321x x x x x问题描述：这是一个五元线性方程组，需要采用全选主元高斯约当消去法来求解。

程序说明：全选主元高斯约当消去法：AGJDN(A,B,N,M,L,JS)A ——双精度实型二维数组，体积为N ×N ，输入参数。

存放方程组的系数矩阵，返回时将被破坏。

B ——双精度实型二维数组，体积为N ×M ，输入兼输出参数。

调用时存放M 组常数向量；返回M 组解向量。

N ——整型变量，输入参数，方程组阶数。

M ——整型变量，输入参数。

方程组右端常数向量的组数。

L ——整型变量，输出参数。

若返回L=0，说明方程组系数矩阵奇异，求解失败；若L ≠0，表示正常返回。

JS ——整型一维数组，长度为N 。

本子程序的工作数组。

全选主元高斯约当消去法的程序及注解如下：子程序：SUBROUTINE AGJDN(A,B,N,M,L,JS)DIMENSION A(N,N),B(N,M),JS(N)DOUBLE PRECISION A,B,DL=1DO 100 K=1,NQ=0.0DO 10 I=K,NDO 10 J=K,NIF (ABS(A(I,J)).GT.Q) THENQ=ABS(A(I,J))JS(K)=JIS=IEND IF10 CONTINUEIF (Q+1.0.EQ.1.0) THENWRITE(*,20)L=0RETURNEND IF20 FORMA T(1X,'FAIL')DO 30 J=K,ND=A(K,J)A(K,J)=A(IS,J)A(IS,J)=D30 CONTINUEDO 40 J=1,MD=B(K,J)B(K,J)=B(IS,J)B(IS,J)=D40 CONTINUEDO 50 I=1,ND=A(I,K)A(I,K)=A(I,JS(K))A(I,JS(K))=D50 CONTINUEDO 60 J=K+1,N60 A(K,J)=A(K,J)/A(K,K)DO 70 J=1,M70 B(K,J)=B(K,J)/A(K,K)DO 90 I=1,NIF (I.NE.K) THENDO 80 J=K+1,N80 A(I,J)=A(I,J)-A(I,K)*A(K,J)DO 85 J=1,M85 B(I,J)=B(I,J)-A(I,K)*B(K,J)END IF90 CONTINUE100 CONTINUEDO 110 K=N,1,-1DO 110 J=1,MD=B(K,J)B(K,J)=B(JS(K),J)B(JS(K),J)=D110 CONTINUERETURNEND主程序：DIMENSION A(5,5),B(5,1),JS(5)DOUBLE PRECISION A,BDA TAA/5.0,7.0,6.0,5.0,1.0,7.0,10.0,8.0,7.0,2.0,6.0,8.0,10.0,9.0,3.0,5.0,7.0,9.0,10.0,4.0,1.0,2.0,3.0,4.0,5. 0/DA TA B/24.0,34.0,35.0,36.0,15.0/N=5M=1CALL AGJDN(A,B,N,M,L,JS)IF (L.NE.0) THENWRITE(*,10) ((B(I,J),I=1,5),J=1,1)END IF10 FORMAT(1X,4D15.6)END运行结果：问题二：二、 利用追赶法求解如下方程组，并给出详细的程序注解和说明： ⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧-----=⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧∙⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---------------72922206431613715211492316521131655232131165214387654321x x x x x x x x 问题描述：这是一个五条对角线的带型的方程组，问题要求需要用追赶法来解。

方法一：追赶法 程序说明：pendag(a,b,c,d,e,r,u,n)n ——整型变量，输入参数，方程阶数a ——n 个元素的一维实型数组，输入参数，存放系数矩阵的下次对角元素下面的元素（,,,321a a a ...n a ),其中21,a a 任意b ——n 个元素的一维实型数组，输入参数，存放系数矩阵的下次对角元素（,,,321b b b ...n a ),其中1b 任意c ——n 个元素的一维实型数组，输入参数，存放系数矩阵的对角元素d ——n 个元素的一维实型数组，输入参数，存放系数矩阵的上次对角元素（,,,321d d d ...n d ),其中n d 任意e ——n 个元素的一维实型数组，输入参数，存放系数矩阵的上次对角元素上面的元素（,,21e e ...,n n n e e e ,,12--)，其中n n e e ,1-任意r ——n 个元素的一维实型数组，输入参数，存放方程的右端向量 u ——n 个元素的一维实型数组，输出参数，输出方程的解向量追赶法带型的计算程序及注解：SUBROUTINE pendag(a,b,c,d,e,r,u,n) PARAMETER(nmax=100)REAL a(n),b(n),c(n),d(n),e(n),r(n),u(n),w(nmax),beta(nmax),alpha(nmax),cg(nmax),h(nmax) INTEGER k,n w(1)=c(1) Beta(1)=0.0Beta(2)=d(1)/w(1) alpha(1)=0.0alpha(2)=e(1)/w(1) alpha(n)=0.0 alpha(n+1)=0.0 Do k=2,nCg(k)=b(k)-a(k)*beta(k-1)W(k)=c(k)-a(k)*alpha(k-1)-cg(k)*beta(k) If(w(k).eq.0.0) pause 'w(k)=0.0 in pendag' Beta(k+1)=(d(k)-cg(k)*alpha(k))/w(k)alpha(k+1)=e(k)/w(k)End doH(1)=0.0H(2)=r(1)/w(1)Do k=2,nH(k+1)=(r(k)-a(k)*h(k-1)-cg(k)*h(k))/w(k)End doU(n)=h(n+1)U(n-1)=h(n)-beta(n)*u(n)Do k=n-2,1,-1U(k)=h(k+1)-beta(k+1)*u(k+1)-alpha(k+1)*u(k+2)End doEnd SUBROUTINE pendagPROGRAM D1R4! Driver program for routine PENDAGPARAMETER(N=8)DIMENSION A1(N,N),A(N),B(N),C(N),D(N),E(N),R(N),U(N),X(N) DA TA A1/3.,-2.,1.,0.,0.,0.,0.,0.,-4.,-5.,3.,2.,&0.,0.,0.,0.,1.,6.,-1.,5.,-3.,0.,0.,0.,0.,1.,2.,-5.,1.,6.,0.,0.,&0.,0.,-3.,6.,-1.,1.,-4.,0.,0.,0.,0.,-1.,2.,-3.,1.,5.,0.,0.,0.,0.,&-5.,2.,-1.,1.,0.,0.,0.,0.,0.,-9.,2.,-7./DA TA R/13.,-6.,-31.,64.,-20.,-22.,-29.,7./Print*,'已知的方程的右端向量'Do I=1,NWRITE(*,'(1X,3F12.6)') R(I)END DODO I=3,NA(I)=A1(I,I-2)END DODO I=2,NB(I)=A1(I,I-1)END DODO I=1,N-1D(I)=A1(I,I+1)END DODO I=1,N-2E(I)=A1(I,I+2)END DODO I=1,NC(I)=A1(I,I)END DOCall PENDAG(A,B,C,D,E,R,U,N)WRITE(*,*)Print*,'计算出的方程的解'DO I=1,NWRITE(*,'(1X,3F12.6)')U(I)END DO!将计算的解B乘以系数矩阵，以检验计算结果的正确DO L=1,NX(L)=0.DO J=1,NX(L)=X(L)+A1(L,J)*U(J)END DOEND DOWRITE(*,*)Print*,'计算出的解乘以系数矩阵的结果'DO I=1,NWRITE(*,'(1X,3F12.6)') X(I)END DOEND运行结果：总结与思考：用完了追赶法，我又偶然发现这个也可以用第一题类似的算法——列选主元高斯消去法算带型问题(方法二）：方法二：列选主元高斯消去法算带型问题程序说明：ABAND(B,D,N,L,IL,M,IT)B——双精度实型二维数据，体积为N×B，输入参数。