fname=程序设计实践.doc

其中的aij是变量xj在第i行的系数，b1,b2,bm等等是已知的常数，而x1，x2，xn等等则是要求的未知数。用线性代数来描述的话，则线性方程组可以写成：AX=B这里的A 是mn 矩阵，为方程组的系数矩阵，X是由n未知数构成的列向量， 是列向量，B 是含有m 个常数列向量。 现在要解决的问题是在已知矩阵 A mn和向量 B m1的情况求得未知向量 X的值。如果有一组数x1、x2、xn使得方程组两边的等号都成立，那么这组数就叫做方程组的解。一个线性方程组的所有的解的集合会被简称为解集。根据解的存在情况，线性方程组可以分为三类: 有唯一解的恰定方程组。 解不存在的超定方程组。 有无穷多解的欠定方程组（也被通俗地称为不定方程组）。下面按照线性代数的形式，分别用前面所说的三种方法，用FORTRAN语言来实现线性方程的求解。5.2.1 高斯消去法高斯消去法实际上就是我们俗称的加减消元法，因著名的数学家高斯和约当得名，它是线性代数中的一个算法，用于决定线性方程组的解，决定矩阵的秩，以及决定可逆方矩阵的逆。当用于一个矩阵时，高斯消去产生“行消去梯形形式”，是一种比较古老的方法，大家在初中数学学习中就接触到了。1.算法原理首先，将前面所述的线性方程组的系数矩阵Amn和常数矩阵B拼成一个增广矩阵，不妨记为Z，形式如下：为了描述方便，在每个变量上加一个上标（i），代表高斯消去法的第i步矩阵中各个数据的值。原始方程组中i=0，因此，增广矩阵Z改写为：假设mn，高斯消去法计算过程如下：（1）从a11开始，通过矩阵的初等变换，将第一列的其他元素变为0；（2）再从新的a22开始重复（1）的过程；（3）以此类推，最多通过m步，最终将增广矩阵变为如下形式：其中，主对角线左面的值都为0。（4）从amm开始，到a22利用矩阵的初等变换，将增广矩阵中第m列到第二列其他元素变为0；（5）每行除以主对角线上的元素，即将对角线上的元素变为1；（6）B列数据即为方程组的解。当m=n且方程组有唯一解时，高斯消去法能很好的求出方程的解。2.程序代码(1)高斯消去法模块：其中：A(N,N)系数矩阵b(N)右向量N方程维数x 方程的根module gausscontainssubroutine solve(A,b,x,N)implicit real*8(a-z)integer:i,k,Nreal*8:A(N,N),b(N),x(N)real*8:Aup(N,N),bup(N)!Ab为增广矩阵 Abreal*8:Ab(N,N+1)Ab(1:N,1:N)=AAb(:,N+1)=b! 这段是 高斯消去法的核心部分do k=1,N-1 do i=k+1,N temp=Ab(i,k)/Ab(k,k) Ab(i,:)=Ab(i,:)-temp*Ab(k,:) end doend do! 经过上一步，增广矩阵已经化为每行主对角线左面的数据都为0的形式。Aup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)!调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solvesubroutine uptri(A,b,x,N)implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1 x(i)=b(i) do j=i+1,N x(i)=x(i)-a(i,j)*x(j) end do x(i)=x(i)/A(i,i)end doend subroutine uptriend module gauss!主程序program mainuse gaussimplicit real*8(a-z)integer,parameter: N=4integer:i,jreal*8:A(N,N),b(N),x(N)!系数矩阵数据存放在h:fsfin.txt中。open(unit=11,file=h:fsfin.txt)!结果存放在h:jgfout.txt中open(unit=12,file=h:jgfout.txt)read(11,*)!读入A矩阵read(11,*)(A(i,j),j=1,N),i=1,N)!读入B向量read(11,*) bcall solve(A,b,x,N)write(12,101)x101 format(T5,高斯消去法计算结果,/,T4,x=,4(/F12.8)print*,b !在屏幕上显示方程的跟。end program main3.运算实例5.2.2 选主元消去法1.算法原理上小节所述的高斯消去法比较简单，但在计算过程中，如果出现下列两种情况：（1）主元为0，算法就没办法继续下去;（2）主元绝对值很小，如小于10-15，以该主元为基础消去，可能使方程组出现很大的误差，方程解变得不太稳定。其中，所谓主元，是指在第k步消去过程中，主对角线上的元素akk。如何解决上述问题，理论上讲，可以选用增广矩阵中主元所在列中非0，的且绝对值较大的任意元素（最后一列的元素除外）来代替主元，比较简单的方式就是选定绝对值最大的元素来代替主元。算法描述如下：（1）从第一列开始，找绝对值最大的元素所在元素，假设为amax11，将第max1行与第一行交换，得到新的a11，通过矩阵的初等变换，将第一列的其他元素变为0；（2）对第二列找除除第一行外，找绝对值最大的元素所在元素，假设为amax21，重复（1）的过程；（3）以此类推，最多通过m步，最终将增广矩阵变为如下形式：其中，主对角线左面的值都为0。（4）从amm开始，到a22利用矩阵的初等变换，将增广矩阵中第m列到第二列其他元素变为0；（5）每行除以主对角线上的元素，即将对角线上的元素变为1；（6）B列数据即为方程组的解。在计算过程中如果出现如下特殊情况，相应处理如下：I、如果某列假设为第k列所有元素都为0，则说明变量xk取任意值，方程都成立，因此方程将有无穷多组解。II、如果该行除最后一列元素外，其余元素都为0，则说明该线性方程组无解，算法结束。2.程序代码module m_gausscontainssubroutine solve(A,b,x,N)implicit real*8(a-z)integer:i,k,Ninteger:id_max !主元素标号real*8:A(N,N),b(N),x(N)real*8:Aup(N,N),bup(N)!Ab为增广矩阵 Abreal*8:Ab(N,N+1)real*8:vtemp1(N+1),vtemp2(N+1)Ab(1:N,1:N)=AAb(:,N+1)=b! 这段是 列主元消去法的核心部分do k=1,N-1 elmax=dabs(Ab(k,k) id_max=k !这段为查找主元素 !这段程序的主要目的不是为了赋值最大元素给elmax，而是为了找出最大元素对应的标号do i=k+1,n if (dabs(Ab(i,k)elmax) then elmax=Ab(i,k) id_max=i end if end do !至此，已经完成查找最大元素，查找完成以后与 第k行交换 !交换两行元素，其他不变 vtemp1=Ab(k,:) vtemp2=Ab(id_max,:) Ab(k,:)=vtemp2 Ab(id_max,:)=vtemp1 do i=k+1,N temp=Ab(i,k)/Ab(k,k) Ab(i,:)=Ab(i,:)-temp*Ab(k,:) end doend doAup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)!调用用上三角方程组的回带方法call uptri(Aup,bup,x,n)end subroutine solvesubroutine uptri(A,b,x,N)implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1 x(i)=b(i) do j=i+1,N x(i)=x(i)-a(i,j)*x(j) end do x(i)=x(i)/A(i,i)end doend subroutine uptriend module m_gaussprogram mainuse m_gaussimplicit real*8(a-z)integer,parameter: N=6integer:i,jreal*8:A(N,N),b(N),x(N)open(unit=11,file=fin.txt)open(unit=12,file=fout.txt)read(11,*)!读入A矩阵read(11,*)(A(i,j),j=1,N),i=1,N)!读入B向量read(11,*) bcall solve(A,b,x,N)write(12,101)x101 format(T5,高斯列主元消去法计算结果,/,T4,x=,6(/F12.8)end program mainimplicit real*8(a-z)integer:N,i,jreal*8:A(N,N),invA(N,N)read(11,*)(A(i,j),j=1,N),i=1,N)write(12,101)101 format(/,T12,逆矩阵为,/)call solve(A,invA,N)write(12,102)(invA(i,j),j=1,n),i=1,n)!变量输出格式只针对 IVF编译器，在CVF中不支持102 format(F16.10/)end subroutine drivesubroutine solve(A,invA,N)implicit real*8(a-z)integer:ninteger:ireal*8:A(n,n),invA(n,n),E(n,n)E=0!设置E为单位矩阵do i=1,n E(i,i)=1end docall mateq(A,E,invA,N,N)end subroutine solvesubroutine mateq(A,B,X,N,M)implicit real*8(a-z)integer:i,k,N,Minteger:id_max !主元素标号real*8:A(N,N),B(N,M),X(N,M)real*8:Aup(N,N),Bup(N,M)!Ab为增广矩阵 ABreal*8:AB(N,N+M)real*8:vtemp1(N+M),vtemp2(N+M)real*8:vtmp(N),xtmp(N)AB(1:N,1:N)=AAB(:,N+1:N+M)=Bdo k=1,N-1 elmax=dabs(Ab(k,k) id_max=kdo i=k+1,n if (dabs(Ab(i,k)elmax) then elmax=Ab(i,k) id_max=i end if end do vtemp1=Ab(k,:) vtemp2=Ab(id_max,:) Ab(k,:)=vtemp2 Ab(id_max,:)=vtemp1 do i=k+1,N temp=Ab(i,k)/Ab(k,k) Ab(i,:)=Ab(i,:)-temp*Ab(k,:) end doend doAup(:,:)=AB(1:N,1:N)do i=1,mvtmp=AB(:,N+i)call uptri(Aup,vtmp,xtmp,n)X(:,i)=xtmpend doend subroutine mateqsubroutine uptri(A,b,x,N)implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)do i=n-1,1,-1 x(i)=b(i) do j=i+1,N x(i)=x(i)-a(i,j)*x(j) end do x(i)=x(i)/A(i,i)end doend subroutine uptriend module inv_matprogram mainuse inv_matinteger:N!N表示矩阵的维数，由文件读入open(unit=11,file=fin.txt)open(unit=12,file=fout.txt)read(11,*)!读说明文字!读入方程维数系数read(11,*)Ncall drive(N)end program main3.运算实例1-范数：；2-范数（欧氏范数）：，1是AHA的最大特征值（AH是A的共轭矩阵，即将A转置，并将A中的虚部取反，如果A是实数矩阵，则AH就是A的转置矩阵AT或记为A）；范数：。范数的一个重要应用就是矩阵的条件数：线性方程组AX=B，其中A为n阶可逆方阵，当B有微小误差B时，考虑解X的误差X。由方程就有：A（X+X）=B+B得，AX=B，因此X=A-1.B，根据范数的定义有：|X|=| A-1.B | A-1|.|B|（1）（注意：这里的推导过程请参看高等工程数学相关书籍的有关章节，下同）。又据AX=B和范数的定义知：|B|A|.|X|，有：（2）将（1）和（2）合并：（3）在（3）式中，左面为X的相对误差，称为方阵A的条件数，记作：Cond（A）条件数事实上表示了矩阵计算对于误差的敏感性。对于线性方程组AX=B，如果A的条件数大，B的微小改变就能引起解X较大的改变，数值稳定性差。如果A的条件数小，B有微小的改变，X的改变也很微小，数值稳定性好。它也可以表示B不变，而A有微小改变时，x的变化情况。比如线性方程组1 2 x = 4 3.999 1 y 7.999的解是(x,y)=(2,1),而1 2 x = 4.001 3.999 1 y 7.998的解是(x,y)=(-3.999,4.000)可见b很小的扰动就引起了x很大的变化，这就是A矩阵条件数大的表现。一个极端的例子，当A奇异时，条件数为无穷，这时即使不改变b，x也可以改变。奇异的本质原因在于矩阵有0特征值，x在对应特征向量的方向上运动不改变Ax的值。如果一个特征值比其它特征值在数量级上小很多，x在对应特征向量方向上很大的移动才能产生b微小的变化，这就解释了为什么这个矩阵为什么会有大的条件数，事实上，正规阵在二范数下的条件数就可以表示成 abs(最大特征值/最小特征值)。我们称B的很小扰动就引起X很大变化的线性方程组AX=B为病态方程。下面给出的是实数矩阵A的2-范数（欧氏范数）和以此为基础的条件数计算的FORTRAN程序设计过程。1.算法原理（1）求A的共轭矩阵AH，因为A为实数矩阵AH就为A的转置矩阵AT。（2）计算AHA。（3）求出AHA的特征值，找出其中最大者并开方即得到2-范数（欧氏范数）：。（4）对于可逆方阵A，求出A的逆矩阵A-1。（5）重复（1）（3）求出A-1的2-范数| A-1|2。（6）计算可逆方阵A的的条件数Cond（A）=|A|2.| A-1|2。主特征值的程序：containssubroutine solve(A,N,namda,u,tol)implicit real*8(a-z)integer:n,i,kreal*8:A(n,n)real*8:u(n),u0(n),v(n)do i=1,n u0(n)=1d0end dou=u0m0=0do k=1,500v=matmul(A,u)call max_rou(v,n,m1)u=v/m1if (dabs(m1-m0)ma) then ma=dabs(r(i) k=i !用k 记录指标，但是ma 不取r(i) end ifend doma=r(k)end subroutine max_rouend module powerprogram mainuse powerimplicit real*8(a-z)real*8:A(3,3),u(3)open(unit=11,file=result.txt)a=reshape(/-1d0,2d0,1d0,& 2d0,-4d0,1d0,& 1d0,1d0,-6d0/),(/3,3/)call solve(A,3,namda,u,1d-7)write(11,101)namda,u101 format(/,T4,幂法计算结果,/,& T3,主特征值为：,/,F12.7,/,& T3,主特征向量为：,3(/F12.7)end program main5.5数据插值与拟合在数学上插值法是函数逼近的一种重要方法，是数值计算的基本课题。而在科学研究中，经常要用到插值法，如：在研究气体的质量M和体积V一定的情况下，压力P与与温度T的关系时，通过实验我们只能得到一组离散的数据，如何通过这些离散数据推断新的数据，插值法就能很好的解决这个问题。插值法根据函数自变量的个数分一维、二维和多维插值等，如上述气体压力的研究中P做因变量，T做自变量，M和V为常数，通过实验方法得到的离散数据，来推断新数据，就是一维插值法。如果除T外，M也当作一个自变量，则用的就是二维插值法。就一维插值法为例，方法比较多，有：多项式插值、lagrange插值、newton插值、hermite插值、分段及样条插值等。在实验数据的处理分析时，数据拟合是经常采用的一种方法，数据拟合的目的是找出能反映变量之间关系的一种表达式，使其在某种准则下最佳的接近已知数据。其原理有最小二乘法、契比雪夫法等，且以最小二乘法最为常见和重要。随着计算技术的快速发展，应用曲线拟合的方法揭示变量之间的联系具有重要的理论和现实意义。本节只讨论一维多项式插值法和最小二乘法拟合以及如何用FORTRAN进行编程计算问题。最小二乘拟合法1.算法原理最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。2.程序代码module lsqcurvefitcontains subroutine solve(x,y,N,c,m)implicit real*8(a-z)integer:N,mreal*8:x(n),y(n),c(M)real*8:bv(M)integer:ireal*8:A(N,M)if (mn) thenwrite(11,101)stopend if101 format(/,warning:the order of polynomial+1 the ,/,& number of date,that is forbidden)do i=1,n call basefunc(bv,x(i),m) A(i,:)=bvend docall leas_eq(A,y,c,N,m)end subroutine solvesubroutine basefunc(bv,x,m)implicit real*8(a-z)integer:Mreal*8:bv(m)integer:ibv(1)=1d0do i=2,m bv(i)=bv(i-1)*xend doend subroutine basefuncsubroutine leas_eq(A,b,x,M,N)implicit real*8(a-z)integer:M,Nreal*8:A(M,N),Q(M,N),R(N,N)real*8:b(M)real*8:QT(N,M) !Q的转置矩阵real*8:QTb(N) !Qbreal*8:x(N)call gram_dec(A,Q,R,M,N)QT=transpose(Q)QTb=matmul(QT,b) ! Rx=Qbcall uptri(R,QTb,x,N) !回带end subroutine leas_eqsubroutine gram_dec(A,Q,R,M,N)implicit real*8(a-z)integer:M,Ninteger:i,j,kreal*8:A(M,N),Q(M,N),R(N,N