最小二乘法应用

最小二乘法应用第一页，共八十五页，编辑于2023年，星期三主要内容1.基本原理2.程序编写3.应用举例第二页，共八十五页，编辑于2023年，星期三AX=b一般m*n线性方程组(m>n)1.基本原理第三页，共八十五页，编辑于2023年，星期三超定方程组：方程的个数超过未知量的个数，一般来说是没有解的。即对于任意一组数(x1,x2,…,xn)不会全为零。我们设想去求一组数第四页，共八十五页，编辑于2023年，星期三使得：取得最小值利用多元函数求极值的方法可得下式第五页，共八十五页，编辑于2023年，星期三即：用矩阵形式给出，即第六页，共八十五页，编辑于2023年，星期三可以证明如果如果A是列满秩的，则该方程存在唯一解，且使得取得最小值，该解就称为超定方程组的最小二乘解。第七页，共八十五页，编辑于2023年，星期三例1：用最小二乘法求下列超定

方程组的近似解第八页，共八十五页，编辑于2023年，星期三解：第九页，共八十五页，编辑于2023年，星期三故得方程组解得x1=0.79272,x2=-1.4641第十页，共八十五页，编辑于2023年，星期三2.程序编写programleast_squareuseIMSLinteger,parameter::m=5integer,parameter::n=2integer::idoubleprecision::A(m,n),B(m),X(n)doubleprecision::C(n,m),D(n,n),E(n)

dataA/2,8,2,7,4,-1,4,1,-1,0/dataB/1,0,1,8,3/

第十一页，共八十五页，编辑于2023年，星期三

if(m>n)thenC=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*Eelseif(m==n)thenX=A.ix.Belsewrite(*,*)"Noroots!"endif

第十二页，共八十五页，编辑于2023年，星期三

doi=1,nwrite(*,*)X(i)enddo

stopendprogramleast_square第十三页，共八十五页，编辑于2023年，星期三如何拟合线性y=ax+b如：有6个数据点，(1,6.9),(2,9.1),(3,10.8)(4,13.2),(5,14.9),(6,17.3)怎么编程求参数a和b？第十四页，共八十五页，编辑于2023年，星期三programleast_squareuseIMSLinteger,parameter::m=6integer,parameter::n=2integer::idoubleprecision::A(m,n),B(m),X(n)doubleprecision::C(n,m),D(n,n),E(n)

dataA/1,2,3,4,5,6,1,1,1,1,1,1/dataB/6.9,9.1,10.8,13.2,14.9,17.3/第十五页，共八十五页，编辑于2023年，星期三if(m>n)thenC=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*Eelseif(m==n)thenX=A.ix.Belsewrite(*,*)"Noroots!"endif第十六页，共八十五页，编辑于2023年，星期三

doi=1,nwrite(*,*)X(i)enddo

stopendprogramleast_square第十七页，共八十五页，编辑于2023年，星期三Ps.in3.应用举例:O2饱和压的拟合与计算55 0.001785755.5 0.002080256 0.00241656.5 0.00279857 0.003231457.5 0.003721758 0.00427558.5 0.00489859 0.005597759.5 0.006381760 0.007258260.5 0.00823661 0.009324361.5 0.01053362 0.01187362.5 0.01335663 0.01499363.5 0.01679764 0.01878164.5 0.020961第一列为温度，单位为K；第二列为实验的饱和压，单位为bar.第十八页，共八十五页，编辑于2023年，星期三65 0.02334965.5 0.02596466 0.02881966.5 0.03193467 0.03532667.5 0.03901568 0.04301968.5 0.0473669 0.05205969.5 0.05713970 0.06262370.5 0.06853571 0.07490171.5 0.08174672 0.08909872.5 0.09698473 0.1054373.5 0.1144874 0.1241474.5 0.13446第十九页，共八十五页，编辑于2023年，星期三75 0.1454775.5 0.1572176 0.1696976.5 0.1829777 0.1970877.5 0.2120578 0.2279278.5 0.2447479 0.2625379.5 0.2813580 0.3012380.5 0.3222281 0.3443681.5 0.3676982 0.3922682.5 0.4181283 0.445383.5 0.4738684 0.5038584.5 0.53532第二十页，共八十五页，编辑于2023年，星期三85 0.5683185.5 0.6028786 0.6390686.5 0.6769387 0.7165387.5 0.7579288 0.8011488.5 0.8462589 0.8933189.5 0.9423890 0.993590.5 1.046791 1.102291.5 1.159892 1.219792.5 1.28293 1.346793.5 1.413994 1.483694.5 1.5559第二十一页，共八十五页，编辑于2023年，星期三95 1.630895.5 1.708596 1.788996.5 1.872297 1.958497.5 2.047598 2.139798.5 2.23599 2.333499.5 2.4351100 2.54100.5 2.6484101 2.7601101.5 2.8753102 2.9941102.5 3.1165103 3.2426103.5 3.3725104 3.5062104.5 3.6438第二十二页，共八十五页，编辑于2023年，星期三105 3.7853105.5 3.9309106 4.0806106.5 4.2344107 4.3925107.5 4.5549108 4.7217108.5 4.8929109 5.0687109.5 5.249110 5.434110.5 5.6237111 5.8183111.5 6.0176112 6.222112.5 6.4313113 6.6458113.5 6.8654114 7.0902114.5 7.3204第二十三页，共八十五页，编辑于2023年，星期三115 7.5559115.5 7.7968116 8.0433116.5 8.2954117 8.5532117.5 8.8166118 9.0859118.5 9.3611119 9.6423119.5 9.9295120 10.223120.5 10.522121 10.828121.5 11.14122 11.459122.5 11.784123 12.115123.5 12.453124 12.798124.5 13.15第二十四页，共八十五页，编辑于2023年，星期三125 13.509125.5 13.874126 14.247126.5 14.627127 15.014127.5 15.408128 15.809128.5 16.218129 16.635129.5 17.059130 17.491130.5 17.93131 18.378131.5 18.833132 19.296132.5 19.768133 20.248133.5 20.736134 21.232134.5 21.737第二十五页，共八十五页，编辑于2023年，星期三135 22.25135.5 22.773136 23.303136.5 23.843137 24.392137.5 24.95138 25.516138.5 26.093139 26.678139.5 27.273140 27.878140.5 28.492141 29.116141.5 29.75142 30.394142.5 31.049143 31.713143.5 32.388144 33.074144.5 33.77第二十六页，共八十五页，编辑于2023年，星期三145 34.477145.5 35.196146 35.925146.5 36.666147 37.418147.5 38.182148 38.958148.5 39.746149 40.547149.5 41.36150 42.186150.5 43.025151 43.878151.5 44.745152 45.626152.5 46.522153 47.434153.5 48.362154 49.307154.5 50.271第二十七页，共八十五页，编辑于2023年，星期三要求：实验数据共200个，要求从文件（文件名：Ps.in)读入，拟合参数，计算出饱和压与实验数据间的误差，找出最大误差和平均误差，并以文件（文件名：Ps.out）形式输出。温度、饱和压及相应中间变量都采取双精度实数。第二十八页，共八十五页，编辑于2023年，星期三氧气饱合压拟合公式如下：其中：c(1)=？c(2)=？编程计算c(3)=？c(4)=？

Tc=154.581!unit:kPc=50.43!unit:bar第二十九页，共八十五页，编辑于2023年，星期三输出文件Ps.out中部分结果

T(K)Psexp(bar)PscalPserr%55.000.00180.0018-0.15755.500.00210.0021-0.12956.000.00240.0024-0.09956.500.00280.0028-0.07457.000.00320.0032-0.05457.500.00370.0037-0.03458.000.00430.0043-0.01658.500.00490.0049-0.00159.000.00560.00560.01159.500.00640.00640.02260.000.00730.00730.03160.500.00820.00820.03861.000.00930.00930.04461.500.01050.01050.05162.000.01190.01190.05462.500.01340.01340.05363.000.01500.01500.05463.500.01680.01680.05664.000.01880.01880.05964.500.02100.02100.054第三十页，共八十五页，编辑于2023年，星期三

150.0042.186042.1836-0.006150.5043.025043.0229-0.005151.0043.878043.8757-0.005151.5044.745044.7424-0.006152.0045.626045.6234-0.006152.5046.522046.5193-0.006153.0047.434047.4308-0.007153.5048.362048.3586-0.007154.0049.307049.3043-0.005154.5050.271050.27080.000

Pserr_ave%=0.015Pserr_max%=-0.157第三十一页，共八十五页，编辑于2023年，星期三程序!ThisprogramistosimulatesaturationpressureofO2.ModuleParametersinteger,parameter::n=200EndModuleParametersProgrammainuseparametersuseIMSLimplicitnonedoubleprecision::T(n),Ps(n),Pscal(n),lnPs(n),thelta(n)doubleprecision::Tc,Pcdoubleprecision::weight(n),coe(4,4),coe1(4),c(4)doubleprecision::Pserr_ave,Pserr_maxinteger::idoubleprecision::Tx,theltax,Pscalx

第三十二页，共八十五页，编辑于2023年，星期三

Tc=154.581!unit:kPc=50.43!unit:baropen(unit=10,file='Ps.in')doi=1,n read(10,*)T(i),Ps(i) weight(i)=1thelta(i)=1-T(i)/Tc lnPs(i)=dlog(Ps(i)/Pc)*T(i)/Tcenddo

callxishu(thelta,lnPs,coe,coe1,weight)

CALLAGAUS(COE,COE1,4,c)!c=coe.ix.coe1

open(unit=20,file='parameters.out')doi=1,4write(20,*)c(i)enddoclose(20)第三十三页，共八十五页，编辑于2023年，星期三

open(unit=30,file='Ps.out')write(30,'(4a12)')"T(K)","Psexp(bar)","Pscal","Pserr%"Pserr_ave=0.Pserr_max=0.

doi=1,n

Pscal(i)=Tc/T(i)*(c(1)*thelta(i)+c(2)*thelta(i)**(1.5)&+c(3)*thelta(i)**(2.5)+c(4)*thelta(i)**5)Pscal(i)=Pc*dexp(Pscal(i))Pserr_ave=Pserr_ave+dabs((Pscal(i)-Ps(i))/Ps(i)*100)if(dabs((Pscal(i)-Ps(i))/Ps(i)*100)>dabs(Pserr_max))thenPserr_max=(Pscal(i)-Ps(i))/Ps(i)*100endif

write(30,'(f12.2,2f12.4,f12.3)')T(i),Ps(i),Pscal(i),(Pscal(i)-Ps(i))/Ps(i)*100

enddo第三十四页，共八十五页，编辑于2023年，星期三

Pserr_ave=Pserr_ave/nwrite(30,*)write(30,*)write(30,'(a20,f8.3)')"Pserr_ave%=",Pserr_avewrite(30,'(a20,f8.3)')"Pserr_max%=",Pserr_maxwrite(30,*)close(30)第三十五页，共八十五页，编辑于2023年，星期三!------------------tocalculatePsatanyT-------------------------!open(unit=40,file='Ps_T.out')write(40,'(2a12)')"T(K)","Ps(bar)"

doi=1,15Tx=80.+(i-1)*5theltax=1-Tx/TcPscalx=Tc/Tx*(c(1)*theltax+c(2)*theltax**(1.5)+c(3)*theltax**(2.5)+& c(4)*theltax**5)Pscalx=Pc*dexp(Pscalx) write(40,'(f12.2,f12.4)')Tx,Pscalx enddoclose(40)!-----------------------------------------------------------------------------------------stopendprogrammain第三十六页，共八十五页，编辑于2023年，星期三subroutinexishu(thelta,lnPs,coe,coe1,weight)useparametersimplicitnonedoubleprecision::coe(4,4),coe1(4),c(4),b(4,n)doubleprecision::thelta(n),lnPs(n),weight(n)integer::i,k,jdoi=1,nb(1,i)=thelta(i)b(2,i)=thelta(i)**(1.5)b(3,i)=thelta(i)**(2.5)b(4,i)=thelta(i)**5 enddo第三十七页，共八十五页，编辑于2023年，星期三

!COE=0.0!COE1=0.0Dok=1,4Doj=1,4COE(k,j)=0.0enddoCOE1(k)=0.0Enddo第三十八页，共八十五页，编辑于2023年，星期三

Dok=1,4Doj=1,4Doi=1,n COE(k,j)=COE(k,j)+b(j,i)*b(k,i)*weight(i)EnddoEnddoDoi=1,n COE1(k)=COE1(k)+lnPs(i)*b(k,i)*weight(i)EnddoEnddo!coe=matmul(b,transpose(b))!coe1=matmul(b,lnps)returnendsubroutinexishu第三十九页，共八十五页，编辑于2023年，星期三!SubroutinetosolvealinearequationgroupSubroutineAGAUS(A,B,N,X)DIMENSIONA(N,N),X(N),B(N),JS(N)DOUBLEPRECISIONA,B,X,T,D

L=1DO50K=1,N-1D=0.0DO210I=K,NDO210J=K,NIF(dABS(A(I,J)).GT.D)THEN D=dABS(A(I,J)) JS(K)=J IS=I ENDIF第四十页，共八十五页，编辑于2023年，星期三210CONTINUEIF(D<=1d-80)THEN L=0 ELSE IF(JS(K).NE.K)THEN DO220I=1,N T=A(I,K) A(I,K)=A(I,JS(K)) A(I,JS(K))=T220CONTINUEENDIF IF(IS.NE.K)THEN DO230J=K,N T=A(K,J) A(K,J)=A(IS,J) A(IS,J)=T第四十一页，共八十五页，编辑于2023年，星期三230CONTINUET=B(K) B(K)=B(IS) B(IS)=T ENDIF ENDIF IF(L.EQ.0)THEN WRITE(*,100) RETURN ENDIF DO10J=K+1,N A(K,J)=A(K,J)/A(K,K)10CONTINUEB(K)=B(K)/A(K,K) DO30I=K+1,N DO20J=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J)第四十二页，共八十五页，编辑于2023年，星期三20CONTINUEB(I)=B(I)-A(I,K)*B(K)30CONTINUE50CONTINUEIF(dABS(A(N,N))<=1d-80)THEN L=0 WRITE(*,100) RETURN ENDIF X(N)=B(N)/A(N,N) DO70I=N-1,1,-1 T=0.0 DO60J=I+1,N T=T+A(I,J)*X(J)60CONTINUE第四十三页，共八十五页，编辑于2023年，星期三

X(I)=B(I)-T70CONTINUE100FORMAT(1X,'FAIL')JS(N)=N DO150K=N,1,-1 IF(JS(K).NE.K)THEN T=X(K) X(K)=X(JS(K)) X(JS(K))=T ENDIF150CONTINUERETURNENDSubroutineAGAUS第四十四页，共八十五页，编辑于2023年，星期三作业拟合y=ax+b直线中的参数a和b，数据点如下：

x,y1，6.92,9.13,10.84,13.25,14.96,17.37,19.18,20.8第四十五页，共八十五页，编辑于2023年，星期三要求数据点以文件形式输入；算出每个点误差；结果也以文件的形式输出。以上面的程序为参考第四十六页，共八十五页，编辑于2023年，星期三求解非线性方程组如：第四十七页，共八十五页，编辑于2023年，星期三可以调用子程序SubroutineNEQNF(FCN,ERRREL,N,ITMAX,XGUESS,X,FNORM)要求：有几个未知数，等式也有几个。第四十八页，共八十五页，编辑于2023年，星期三例：programmainuseIMSLimplicitnoneexternalFCNreal,parameter::ERRREL=0.0001integer,parameter::N=3integer,parameter::ITMAX=100real::XGUESS(N)=(/0.0,1.0,2.0/)realX(N),FNORMCALLNEQNF(FCN,ERRREL,N,ITMAX,XGUESS,X,FNORM)write(*,*)Xstopend第四十九页，共八十五页，编辑于2023年，星期三subroutineFCN(XA,F,N)implicitnoneintegerNreal,target::XA(N)realF(N)real,pointer::x,y,z!在计算时使用x,y,z看起来比较清楚

x=>XA(1)y=>XA(2)z=>XA(3)F(1)=x*x+y*y+z*z-3F(2)=x*y+y*z+x*z-3F(3)=exp(x)+exp(y)+exp(z)-3*exp(1.0)returnendsubroutine第五十页，共八十五页，编辑于2023年，星期三求解超定非线性方程组？要使用非线性拟合程序，比较复杂。通常有几个未知量，就要给几个初值。往往每给定一组初值，由目标函数得到的最小平均误差也不一样。所以要根据实际情况，进行多次尝试和分析，得出最后的参数值。第五十一页，共八十五页，编辑于2023年，星期三可化为多元线性回归函数式通式：变换变量变换常数yxabcylnxcosxcabxabcx/yxabcy1/xabclnyxlnabc第五十二页，共八十五页，编辑于2023年，星期三作业程序1!Thisprogramistosimulateparametersofliney=ax+bModuleParametersinteger,parameter::n=8EndModuleParametersProgrammainuseparametersuseIMSLimplicitnonedoubleprecision::x(n),y(n),ycal(n)doubleprecision::coe(2,2),coe1(2),c(2)doubleprecision::yerr_ave,yerr_maxinteger::i第五十三页，共八十五页，编辑于2023年，星期三

open(unit=10,file='data.in')doi=1,nread(10,*)x(i),y(i)enddo

callxishu(x,y,coe,coe1)

c=coe.ix.coe1

open(unit=20,file='parameters.out')doi=1,2write(20,*)c(i)enddoclose(20)第五十四页，共八十五页，编辑于2023年，星期三

open(unit=30,file='data.out')write(30,'(4a12)')"x","yexp","ycal","yserr%"yerr_ave=0.yerr_max=0.

doi=1,n

ycal(i)=c(1)*x(i)+c(2)yerr_ave=yerr_ave+dabs((ycal(i)-y(i))/y(i)*100)if(dabs((ycal(i)-y(i))/y(i)*100)>dabs(yerr_max))thenyerr_max=(ycal(i)-y(i))/y(i)*100endif

write(30,'(f12.2,2f12.4,f12.3)')x(i),y(i),ycal(i),(ycal(i)-y(i))/y(i)*100

enddo第五十五页，共八十五页，编辑于2023年，星期三

yerr_ave=yerr_ave/nwrite(30,*)write(30,*)write(30,'(a20,f8.3)')"yerr_ave%=",yerr_avewrite(30,'(a20,f8.3)')"yerr_max%=",yerr_maxwrite(30,*)close(30)stopendprogrammain第五十六页，共八十五页，编辑于2023年，星期三subroutinexishu(x,y,coe,coe1)useparametersimplicitnonedoubleprecision::coe(2,2),coe1(2),c(2),b(2,n)doubleprecision::x(n),y(n)integer::i,k,jdoi=1,nb(1,i)=x(i)b(2,i)=1.0enddo

coe=matmul(b,transpose(b))coe1=matmul(b,y)returnendsubroutinexishu第五十七页，共八十五页，编辑于2023年，星期三作业程序2programleast_squareuseIMSLinteger,parameter::m=8integer,parameter::n=2integer::idoubleprecision::A(m,n),B(m),X(n)doubleprecision::C(n,m),D(n,n),E(n)

dataA/1,2,3,4,5,6,7,8,1,1,1,1,1,1,1,1/dataB/6.9,9.1,10.8,13.2,14.9,17.3,19.1,20.8/第五十八页，共八十五页，编辑于2023年，星期三

if(m>n)thenC=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*Eelseif(m==n)thenX=A.ix.Belsewrite(*,*)"Noroots!"endif第五十九页，共八十五页，编辑于2023年，星期三doi=1,nwrite(*,*)X(i)enddo

stopendprogramleast_square第六十页，共八十五页，编辑于2023年，星期三作业程序3programleast_squareuseIMSLinteger,parameter::m=8integer,parameter::n=2integer::idoubleprecision::A(m,n),B(m),X(n)doubleprecision::C(n,m),D(n,n),E(n)doubleprecision::ycal,yerr_ave,yerr_max第六十一页，共八十五页，编辑于2023年，星期三

open(unit=10,file='data.in')doi=1,mread(10,*)A(i,1),B(i)enddodoi=1,mA(i,2)=1.0enddo第六十二页，共八十五页，编辑于2023年，星期三

if(m>n)thenC=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*Eelseif(m==n)thenX=A.ix.Belsewrite(*,*)"Noroots!"endif

open(unit=20,file='parameters.out')doi=1,nwrite(20,*)X(i)enddo第六十三页，共八十五页，编辑于2023年，星期三

open(unit=30,file='data.out')write(30,'(4a12)')"x","yexp","ycal","yserr%"yerr_ave=0.yerr_max=0.

doi=1,m

ycal=x(1)*A(i,1)+x(2)yerr_ave=yerr_ave+dabs((ycal-B(i))/B(i)*100)if(dabs((ycal-B(i))/B(i)*100)>dabs(yerr_max))thenyerr_max=(ycal-B(i))/B(i)*100endif

write(30,'(f12.2,2f12.4,f12.3)')A(i,1),B(i),ycal,(ycal-B(i))/B(i)*100

enddo第六十四页，共八十五页，编辑于2023年，星期三

yerr_ave=yerr_ave/mwrite(30,*)write(30,*)write(30,'(a20,f8.3)')"yerr_ave%=",yerr_avewrite(30,'(a20,f8.3)')"yerr_max%=",yerr_maxwrite(30,*)close(30)

stopendprogramleast_square第六十五页，共八十五页，编辑于2023年，星期三Rline函数拟合直线y=ax+b直线SubroutineRLINE(NOBS,XDATA,YDATA,B0,B1,STAT)integerNOBS:数据点数目realXDATA(NOBS)：数据点X轴值realyDATA(NOBS)：数据点Y轴值realB0,B1:Y=B0+B1*Xrealstat(12):返回的信息。第六十六页，共八十五页，编辑于2023年，星期三例：programmainuseIMSLimplicitnoneinteger,parameter::NOBS=5realXDATA(NOBS),YDATA(NOBS)realB0,B1realSTAT(12)realxinc,xp,valuerealF,XF(X)=2*X+3integerI第六十七页，共八十五页，编辑于2023年，星期三

!产生数据点

doI=1,NOBSXDATA(I)=real(I)YDATA(I)=F(XDATA(I))enddocallRLINE(NOBS,XDATA,YDATA,B0,B1,STAT)!结果一定为y=2X+3,因为数据点是根据这个函数来产生的.write(*,"('Y=',F5.2,'X+'F5.2)")B1,B0stopendprogram第六十八页，共八十五页，编辑于2023年，星期三Rcurv函数用最小二乘法来求一条最接近数据点的n次多项式。SubroutineRCURV(NOBS,XDATA,YDATA,NDEG,B,SSPOLY,STAT)integerNOBS:数据点数目realXDATA(NOBS)：数据点X轴值realyDATA(NOBS)：数据点Y轴值integerNDEG:所要计算的多项式次数第六十九页，共八十五页，编辑于2023年，星期三RealB(NDEG+1):多项式系数

realSSPOLY:返回信息

realSTAT(10):返回信息第七十页，共八十五页，编辑于2023年，星期三例：programmainuseIMSLimplicitnoneinteger,parameter::NOBS=10,NDEG=2realXDATA(NOBS),YDATA(NOBS)realB(NDEG+1),SSPOLY(NDEG+1)realSTAT(12)realF,X,RF(X,R)=R+1+2*X+3*X*XintegerI第七十一页，共八十五页，编辑于2023年，星期三callrandom_seed()!产生数据点

doI=1,NOBSXDATA(I)=real(I)callrandom_number(R)YDATA(I)=F(XDATA(I),R)enddocallRCURV(NOBS,XDATA,YDATA,NDEG,B,SSPOLY,STAT)write(*,"(F5.2'+'F5.2'X+'F5.2'X*X')")Bstopendprogram第七十二页，共八十五页，编辑于2023年，星期三小结F(x)=0求根

二分法，牛顿迭代法线性方程(n*n)组求解：

高斯消元法，逆矩阵法超定线性方程(m*n,m>n)组求解：

最小二乘法性线和非线性拟合：

最小二乘法第七十三页，共八十五页，编辑于2023年，星期三应用举例如何求范德华方程中气相的体积？第七十四页，共八十五页，编辑于2023年，星期三方程展开以后为：对于水来说：第七十五页，共八十五页，编辑于2023年，星期三应用程序programmainimplicitnonereal::T,P,V!units:K,bar,cm3/molreal::a,breal::V1,V2,V3real::Psreal,parameter::R=83.14472!unit:cm3.bar/(mol.K)real,parameter::Tc=647.096!unit:Kreal,parameter::Pc=220.64!unit:bar第七十六页，共八十五页，编辑于2023年，星期三100write(*,*)write(*,*)"InputT(K)andP(bar):"read(*,*)T,P

if(T<273.16)thenwrite(*,*)"wrongT"goto100endifif(T>=273.16.and.T<=647.096)thencallH2Osaturation(T,Ps)endif

if(T>=273.16.and.T<=647.096.and.P>Ps)thenwrite(*,*)"Notvaporphase!"goto100endif第七十七页，共八十五页，编辑于2023年，星期三

a=27.0*R**2*Tc**2/(64.0*Pc)b=R*Tc/(8.0*Pc)callroots(-(b+R*T/P),a/P,-a*b/P,V1,V2,V3)V=max(V1,V2,V3)write(*,*)write(*,*)"T=","K"write(*,*)"P=","bar"write(*,*)"V=",V,"cm3/mol"write(*,*)stopendprogrammain第七十八页，共八十五页，编辑于2023年，星期三!fromWagnerandPruss,J.Phys.Chem.Ref.Data,22(3),1993.s