Question

关于fortran 90里面的编程问题（急）

设计任务（选择一题）1．线性方程组求解问题。一物理系统可用下列线性方程组来表示：从文件中读入m1、m2和θ的值，求a1、a2、N1 和N2的值。其中g取9.8，输入θ时以角度为单位。要求：（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三角分解法、追赶法等），定义求解线性方程组Ax=b的子程序，要求该子程序能求解任意线性方程组。（2）在主程序中分别调用上面定义的两个子程序，并对求解结果进行对比分析。（3）绘制以上两个方法所求得的方程解的数据分布图。2．线性病态方程组问题。下面是一个线性病态方程组：（1）求方程的解。（2）将方程右边向量元素b3改为0.53，再求解，并比较b3的变化和解的相对变化。（3）计算系数矩阵A的条件数并分析结论。提示：矩阵A的条件数等于A的范数与A的逆矩阵的范数的乘积，即 。这样定义的条件数总是大于1的。条件数越接近于1，矩阵的性能越好，反之，矩阵的性能越差。矩阵A的条件数 ，其中 ，aij系矩阵A的元素。要求：（1）方程的系数矩阵、常数向量均从文件中读入。（2）定义求解线性方程组Ax=b的子程序，要求该子程序能求解任意线性方程组。（3）在主程序中调用子程序，并对求解结果进行对比分析。（4）绘制常数向量修改前后所求得的方程解的数据分布图。3．数值迭代计算问题。已知 ，其中h=0.002，k是迭代次数。 （1）计算 及 ，要求结果精确到小数点后第6位，并将结果输出到文件poisson.dat中。（2）找出 及 的最大、最小值，并统计两个极值各有几个。（3）统计 在0.000002以下、0.000002~0.000029、0.000003~0.000049、0.000005~0.000010、0.000010以上的数据个数，并绘制数据分布饼图或直方图。

Answer 1

subroutine equil (aik,y,fe,elem,nele,ngas,ntot,ind)
c aik is elemental specie composition matrix nele x ntotal dimension
c y is no moles of each specie ntotal dimension
c fe is free energy of each specie ntotal dimension
c elem is no moles of each element nele dimension
c nele is no elements
c ngas is no gas species
c ntot is total no species
c ind is -7 if singular matrix and -1 if other than last solid disappeared
dimension aik(250),y(25),fe(25),elem(10),x(25),bmat(36),amat(1296)
c this common statement special to bkw code
common/subvar/amaxe,aminx,aminy,tx(10)
sst=0.
cpt=0.0
do 2 i=1,ntot
if(y(i).lt.aminy) y(i)=aminy
2 continue
35 nm1=ntot+nele+1
nn1=ntot-ngas
nm1sq=nm1*nm1
c zero amatrix
do 4 i=1,nm1sq
amat(i)=0.
4 continue
c form sums
sum=0.
do 5 i=1,ngas
sum=sum+y(i)
5 continue
bary=sum
rbary=(1.)/bary
c fill b matrix
i=1
do 6 j=1,ngas
bmat(i)=-(fe(i)+alog(rbary*y(i)))
i=i+1
6 continue
if(nn1.eq.0) go to 70
do 7 j=1,nn1
bmat(i)=-(fe(i))
i=i+1
7 continue
70 continue
do 8 j=1,nele
bmat(i)=elem(j)
i=i+1
8 continue
bmat(i)=0.
c fill in amatrix
i=1
do 9 j=1,ngas
amat(i)=(1.)/(y(j))
i=i+nm1+1
9 continue
i=ntot+1
j=1
l=1
12 do 10 k=1,nele
amat(i)=aik(j)
i=i+1
j=j+1
10 continue
if(l.gt.ngas)go to 11
amat(i)=1.0
11 l=l+1
if(l.gt.ntot)go to 13
i=i+1+ntot
go to 12
13 i=i+1
do 14 k=1,ngas
amat(i)=-rbary
i=i+1
14 continue
i=i+nn1+nele
amat(i)=-1.0
i=i+1
l=1
j=1
17 do 15 k=1,ntot
amat(i)=aik(j)
j=j+nele
i=i+1
15 continue
i=i+nele+1
l=l+1
if(l.gt.nele) go to 16
j=l
go to 17
16 call lss(nm1,1,nm1,amat,bmat,d,det,ind)
if(ind.lt.0) go to 40
c put answers in x
do 18 i=1,ntot
x(i)=bmat(i)
18 continue
c test for too small x
do 19 i=1,ngas
if(x(i).lt.aminx) x(i)=aminx
19 continue
c test to see if any solids disappeared
if(nn1.eq.0) go to 22
i=ngas+1
do 20 j=1,nn1
if(x(i).lt.0.) go to 50
i=i+1
20 continue
22 continue
c test to see if converged
80 sum=0.
do 21 i=1,ntot
sum=sum+abs(y(i)-x(i))
21 continue
if(sum.lt.amaxe) go to 60
c special go ahead with less strict error limit
cpt=cpt+1.0
if(cpt.gt.50.)write (3,991)sum
991 format(25h equil amaxe accepted was,1e18.11)
if(cpt.gt.50.)go to 60
c resolve with y now having last answers x
31 do 32 i=1,ntot
y(i)=x(i)
32 continue
go to 35
c have conve-ged
60 do 38 i=1,ntot
y(i)=x(i)
38 continue
j=ntot+1
if(sst.le.0.)go to 34
y(j)=0.
ntot=ntot+1
34 return
c solid has disappeared
50 sst=1.0
if(i.lt.ntot)go to 33
ntot=ntot-1
go to 35
c an error has occurred as other than last solid disappeared
c ind=-1 if error has occurred
33 ind=-1
return
c error has occurred matrix is singular
40 ind=-7
return
end