fortran90解线性方程组程序.docVIP

扬州大学 解线性方程组 计算流体动力学 水研2010级 第 PAGE 5 页 共 NUMPAGES 6 页 水利学院 题目： 用Fortran语言编一程序解下列n元一次方程组，并对相同的线性方程组改变其矢量项的一个参数，对比其计算结果的变化。 解题方法： 采用解线性方程组直接算法中的高斯列主元消去法。 高斯列主元消去法，较高斯消去法舍入误差小，较高斯全主元消去法计算量少，所以这里采用高斯列主元消去法。 算法基本思想： 该方法分消元和回代两个过程，对其增广矩阵进行转换。 1.消元过程：对k=1，2，…,(n-1)，依次计算 其中， 对方程，其增广矩阵，在每一步消元前先进行按列选主元的工作，如第k步消元前，就所有，取绝对值最大值，设，将第L个方程与第 i个方程进行互换位置，这样成为第k步的主元素，然后进行第k步消元。每步消元都如此，最终将其化为上三角矩阵。 2.回代过程： 依次分别计算出Xn，Xn-1，Xn-2，……X2，X1。 源程序如下： program linear_equation_bc !主程序开始语句 implicit none !说明部分 integer::i,j,k,imax,t !定义整型变量，i、j分别标记增广矩阵的行列 !imax标记列主元最大的行号，t为输入数组元数 real::max,n !max记录每次计算的列主元的值，n方程组的元数 !用矩阵(实型动态数组)将线性方程组表示出来 real,dimension(:,:),allocatable::a,m !定义实型二维数组a、m，分别存放初始 !增广矩阵和选主元时交换对应行时的临时矩阵 real,dimension(:),allocatable::x !定义实型一维数组x，存放计算结果 print*,请输入线性方程组的元数n： !表控输出语句 read*,t !输入语句，输入数组的元数 !给动态数组分配内存 allocate(a(t,t+1),m(t,t+1),x(t)) !a、m为t行t+1列二维数组，x为t行一维数组 print*,请依次按行输入线性方程组的增广矩阵a： !表控输出语句 read*,((a(i,j),j=1,t+1),i=1,t) !隐do循环输入语句，输入增广矩阵a do k =1,t-1 !guass列主元消去法消元过程开始 !下面选取每列的最大列主元素 max=abs(a(k,k)) imax=k do i=k+1,t+1 !do语句 if (abs(a(i,k))max) then !do块 !if-then语句 max=abs(a(i,k)) !if-then块 imax=i end if !end if语句 end do !end do语句 !将最大列元素所在的行与第K行进行交换 do j=k,t+1 m(k,j)=a(k,j) a(k,j)=a(imax,j) a(imax,j)=m(k,j) end do !对方程组按X1, X2,…, Xt的顺序进行依次消元 !将增广矩阵a化为上三角矩阵 do i=k+1,t m(i,k)=a(i,k)/a(k,k) do j=k+1,t+1 a(i,j)=a(i,j)-m(i,k)*a(k,j) end do end do end do !guass列主元消去法的第一步消去至此结束 x(t)= a(t,t+1)/a(t,t) !先计算出xt do k=t-1,1,-1 !回代计算开始 n=0 !给定n初值，n作为累加变量，记录每次计算中的累加值 do j=t,k+1,-1 !do循环，计算n的累加值 n=n+a(k,j)*x(j) end do x(k)=( a(k,t+1)-n)/a(k,k) !根据公式算的x(k)的值，即为计算结果 end do !回代过程结束 print*,线性方程组的解为： !表控输出语句 do i =1,t !用显do循环输出计算结果 print*,x(i)