最小二乘法fortran应用

最小二乘法最小二乘法计算地球化学主要内容主要内容n n1.基本原理n n2.程序编写n n3.应用举例AXAX=b b一般一般m*n线性方程组线性方程组(mn)1.基本原理基本原理超定方程组超定方程组：方程的个数超过未知量的个数，一般来说：方程的个数超过未知量的个数，一般来说是没有解的。即对于任意一组数是没有解的。即对于任意一组数(x1,x2,xn)不会全为零。我们设想去求一组数不会全为零。我们设想去求一组数使得：使得：取得最小值取得最小值利用多元函数求极值的方法可得下式利用多元函数求极值的方法可得下式即：即：用矩阵形式给出，即n n可以证明如果如果可以证明如果如果A A是列满秩的，则该方程存在是列满秩的，则该方程存在唯一解，且使得唯一解，且使得取得最小值，该解就称为超定方程组的最小二乘解。例例 1：用最小二乘法求下列超定用最小二乘法求下列超定方程组的近似解方程组的近似解n n解：解：故得方程组故得方程组解得解得 x1=0.79272,x2=-1.46412.程序编写程序编写program least_square use IMSL integer,parameter:m=5 integer,parameter:n=2 integer :i double precision :A(m,n),B(m),X(n)double precision :C(n,m),D(n,n),E(n)data A/2,8,2,7,4,-1,4,1,-1,0 /data B/1,0,1,8,3/if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*E else if(m=n)then X=A.ix.B else write(*,*)No roots!end if do i=1,n write(*,*)X(i)end do stopend program least_square如何拟合线性如何拟合线性y=ax+bn n如：有6个数据点，(1,6.9),(2,9.1),(3,10.8)(4,13.2),(5,14.9),(6,17.3)怎么编程求参数 a 和 b？program least_square use IMSL integer,parameter:m=6 integer,parameter:n=2 integer :i double precision :A(m,n),B(m),X(n)double precision :C(n,m),D(n,n),E(n)data A/1,2,3,4,5,6,1,1,1,1,1,1 /data B/6.9,9.1,10.8,13.2,14.9,17.3/if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*Eelse if(m=n)then X=A.ix.B else write(*,*)No roots!end if do i=1,n write(*,*)X(i)end do stopend program least_squaren nPs.in3.应用举例应用举例:O2饱和压的拟合与计算饱和压的拟合与计算550.001785755.50.0020802560.00241656.50.002798570.003231457.50.0037217580.00427558.50.004898590.005597759.50.0063817600.007258260.50.008236610.009324361.50.010533620.01187362.50.013356630.01499363.50.016797640.01878164.50.020961第一列为温度，单位为第一列为温度，单位为K；第二列为实验的饱和压，单位为第二列为实验的饱和压，单位为bar.650.02334965.50.025964660.02881966.50.031934670.03532667.50.039015680.04301968.50.04736690.05205969.50.057139700.06262370.50.068535710.07490171.50.081746720.08909872.50.096984730.1054373.50.11448740.1241474.50.13446750.1454775.50.15721760.1696976.50.18297770.1970877.50.21205780.2279278.50.24474790.2625379.50.28135800.3012380.50.32222810.3443681.50.36769820.3922682.50.41812830.445383.50.47386840.5038584.50.53532850.5683185.50.60287860.6390686.50.67693870.7165387.50.75792880.8011488.50.84625890.8933189.50.94238900.993590.51.0467911.102291.51.1598921.219792.51.282931.346793.51.4139941.483694.51.5559951.630895.51.7085961.788996.51.8722971.958497.52.0475982.139798.52.235992.333499.52.43511002.54100.52.64841012.7601101.52.87531022.9941102.53.11651033.2426103.53.37251043.5062104.53.64381053.7853105.53.93091064.0806106.54.23441074.3925107.54.55491084.7217108.54.89291095.0687109.55.2491105.434110.55.62371115.8183111.56.01761126.222112.56.43131136.6458113.56.86541147.0902114.57.32041157.5559115.57.79681168.0433116.58.29541178.5532117.58.81661189.0859118.59.36111199.6423119.59.929512010.223120.510.52212110.828121.511.1412211.459122.511.78412312.115123.512.45312412.798124.513.1512513.509125.513.87412614.247126.514.62712715.014127.515.40812815.809128.516.21812916.635129.517.05913017.491130.517.9313118.378131.518.83313219.296132.519.76813320.248133.520.73613421.232134.521.73713522.25135.522.77313623.303136.523.84313724.392137.524.9513825.516138.526.09313926.678139.527.27314027.878140.528.49214129.116141.529.7514230.394142.531.04914331.713143.532.38814433.074144.533.7714534.477145.535.19614635.925146.536.66614737.418147.538.18214838.958148.539.74614940.547149.541.3615042.186150.543.02515143.878151.544.74515245.626152.546.52215347.434153.548.36215449.307154.550.271要求：要求：n n实验数据共200个，要求从文件（文件名：Ps.in)读入，拟合参数，计算出饱和压与实验数据间的误差，找出最大误差和平均误差，并以文件（文件名：Ps.out）形式输出。n n温度、饱和压及相应中间变量都采取双精度实数。氧气饱合压拟合公式如下：氧气饱合压拟合公式如下：其中其中：c(1)=？c(2)=？编程计算c(3)=？c(4)=？Tc=154.581!unit:k Pc=50.43 !unit:bar输出文件输出文件Ps.out中部分结果中部分结果 T(K)Psexp(bar)Pscal Pserr%55.00 0.0018 0.0018 -0.157 55.50 0.0021 0.0021 -0.129 56.00 0.0024 0.0024 -0.099 56.50 0.0028 0.0028 -0.074 57.00 0.0032 0.0032 -0.054 57.50 0.0037 0.0037 -0.034 58.00 0.0043 0.0043 -0.016 58.50 0.0049 0.0049 -0.001 59.00 0.0056 0.0056 0.011 59.50 0.0064 0.0064 0.022 60.00 0.0073 0.0073 0.031 60.50 0.0082 0.0082 0.038 61.00 0.0093 0.0093 0.044 61.50 0.0105 0.0105 0.051 62.00 0.0119 0.0119 0.054 62.50 0.0134 0.0134 0.053 63.00 0.0150 0.0150 0.054 63.50 0.0168 0.0168 0.056 64.00 0.0188 0.0188 0.059 64.50 0.0210 0.0210 0.054 150.00 42.1860 42.1836 -0.006 150.50 43.0250 43.0229 -0.005 151.00 43.8780 43.8757 -0.005 151.50 44.7450 44.7424 -0.006 152.00 45.6260 45.6234 -0.006 152.50 46.5220 46.5193 -0.006 153.00 47.4340 47.4308 -0.007 153.50 48.3620 48.3586 -0.007 154.00 49.3070 49.3043 -0.005 154.50 50.2710 50.2708 0.000 Pserr_ave%=0.015 Pserr_max%=-0.157程程 序序!This program is to simulate saturation pressure of O2.Module Parameters integer,parameter :n=200End Module ParametersProgram mainuse parametersuse IMSLimplicit none double precision :T(n),Ps(n),Pscal(n),lnPs(n),thelta(n)double precision :Tc,Pc double precision :weight(n),coe(4,4),coe1(4),c(4)double precision :Pserr_ave,Pserr_max integer :i double precision :Tx,theltax,Pscalx Tc=154.581!unit:k Pc=50.43 !unit:bar open(unit=10,file=Ps.in)do i=1,n read(10,*)T(i),Ps(i)weight(i)=1 thelta(i)=1-T(i)/TclnPs(i)=dlog(Ps(i)/Pc)*T(i)/Tc end do call xishu(thelta,lnPs,coe,coe1,weight)CALL AGAUS(COE,COE1,4,c)!c=coe.ix.coe1 open(unit=20,file=parameters.out)do i=1,4 write(20,*)c(i)end do close(20)open(unit=30,file=Ps.out)write(30,(4a12)T(K),Psexp(bar),Pscal,Pserr%Pserr_ave=0.Pserr_max=0.do i=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)then Pserr_max=(Pscal(i)-Ps(i)/Ps(i)*100 end if write(30,(f12.2,2f12.4,f12.3)T(i),Ps(i),Pscal(i),(Pscal(i)-Ps(i)/Ps(i)*100 end do Pserr_ave=Pserr_ave/n write(30,*)write(30,*)write(30,(a20,f8.3)Pserr_ave%=,Pserr_ave write(30,(a20,f8.3)Pserr_max%=,Pserr_max write(30,*)close(30)!-to calculate Ps at any T-!open(unit=40,file=Ps_T.out)write(40,(2a12)T(K),Ps(bar)do i=1,15 Tx=80.+(i-1)*5 theltax=1-Tx/Tc Pscalx=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 end do close(40)!-stopend program mainsubroutine xishu(thelta,lnPs,coe,coe1,weight)use parameters implicit none double precision:coe(4,4),coe1(4),c(4),b(4,n)double precision:thelta(n),lnPs(n),weight(n)integer :i,k,j do i=1,n b(1,i)=thelta(i)b(2,i)=thelta(i)*(1.5)b(3,i)=thelta(i)*(2.5)b(4,i)=thelta(i)*5 end do !COE=0.0 !COE1=0.0 Do k=1,4 Do j=1,4 COE(k,j)=0.0 end do COE1(k)=0.0 End do Do k=1,4 Do j=1,4 Do i=1,n COE(k,j)=COE(k,j)+b(j,i)*b(k,i)*weight(i)Enddo Enddo Do i=1,n COE1(k)=COE1(k)+lnPs(i)*b(k,i)*weight(i)Enddo Enddo !coe=matmul(b,transpose(b)!coe1=matmul(b,lnps)returnend subroutine xishu!Subroutine to solve a linear equation group Subroutine AGAUS(A,B,N,X)DIMENSION A(N,N),X(N),B(N),JS(N)DOUBLE PRECISION A,B,X,T,D L=1 DO 50 K=1,N-1 D=0.0 DO 210 I=K,N DO 210 J=K,N IF(dABS(A(I,J).GT.D)THEN D=dABS(A(I,J)JS(K)=J IS=I ENDIF210 CONTINUE IF(D=1d-80)THEN L=0ELSE IF(JS(K).NE.K)THEN DO 220 I=1,N T=A(I,K)A(I,K)=A(I,JS(K)A(I,JS(K)=T220 CONTINUE ENDIF IF(IS.NE.K)THEN DO 230 J=K,N T=A(K,J)A(K,J)=A(IS,J)A(IS,J)=T230 CONTINUE T=B(K)B(K)=B(IS)B(IS)=T ENDIF ENDIF IF(L.EQ.0)THEN WRITE(*,100)RETURN ENDIF DO 10 J=K+1,N A(K,J)=A(K,J)/A(K,K)10 CONTINUE B(K)=B(K)/A(K,K)DO 30 I=K+1,N DO 20 J=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J)20 CONTINUE B(I)=B(I)-A(I,K)*B(K)30 CONTINUE50 CONTINUE IF(dABS(A(N,N)XA(1)y=XA(2)z=XA(3)F(1)=x*x+y*y+z*z-3 F(2)=x*y+y*z+x*z-3 F(3)=exp(x)+exp(y)+exp(z)-3*exp(1.0)returnend subroutine求解超定非线性方程组？求解超定非线性方程组？n n要使用非线性拟合程序，比较复杂。通常有几个未知量，就要给几个初值。往往每给定一组初值，由目标函数得到的最小平均误差也不一样。所以要根据实际情况，进行多次尝试和分析，得出最后的参数值。可化为多元线性回归可化为多元线性回归函数式函数式 通式：通式：变换变量变换变量变换常数变换常数y yx xa ab bc cy ylnxlnxcosxcosxc ca ab bx xa ab bc cx/yx/yx xa ab bc cy y1/x1/xa ab bc clnylnyx xlnalnab bc c作业程序作业程序1!This program is to simulate parameters of line y=ax+bModule Parameters integer,parameter :n=8End Module ParametersProgram mainuse parametersuse IMSLimplicit none double precision :x(n),y(n),ycal(n)double precision :coe(2,2),coe1(2),c(2)double precision :yerr_ave,yerr_max integer :i open(unit=10,file=data.in)do i=1,n read(10,*)x(i),y(i)end do call xishu(x,y,coe,coe1)c=coe.ix.coe1 open(unit=20,file=parameters.out)do i=1,2 write(20,*)c(i)end do close(20)open(unit=30,file=data.out)write(30,(4a12)x,yexp,ycal,yserr%yerr_ave=0.yerr_max=0.do i=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)then yerr_max=(ycal(i)-y(i)/y(i)*100 end if write(30,(f12.2,2f12.4,f12.3)x(i),y(i),ycal(i),(ycal(i)-y(i)/y(i)*100 end do yerr_ave=yerr_ave/n write(30,*)write(30,*)write(30,(a20,f8.3)yerr_ave%=,yerr_ave write(30,(a20,f8.3)yerr_max%=,yerr_max write(30,*)close(30)stopend program mainsubroutine xishu(x,y,coe,coe1)use parametersimplicit none double precision:coe(2,2),coe1(2),c(2),b(2,n)double precision:x(n),y(n)integer :i,k,j do i=1,n b(1,i)=x(i)b(2,i)=1.0 end do coe=matmul(b,transpose(b)coe1=matmul(b,y)returnend subroutine xishu作业程序作业程序2program least_square use IMSL integer,parameter:m=8 integer,parameter:n=2 integer :i double precision :A(m,n),B(m),X(n)double precision :C(n,m),D(n,n),E(n)data A/1,2,3,4,5,6,7,8,1,1,1,1,1,1,1,1 /data B/6.9,9.1,10.8,13.2,14.9,17.3,19.1,20.8/if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*E else if(m=n)then X=A.ix.B else write(*,*)No roots!end ifdo i=1,n write(*,*)X(i)end do stopend program least_square作业程序作业程序3program least_square use IMSL integer,parameter:m=8 integer,parameter:n=2 integer :i double precision :A(m,n),B(m),X(n)double precision :C(n,m),D(n,n),E(n)double precision :ycal,yerr_ave,yerr_max open(unit=10,file=data.in)do i=1,m read(10,*)A(i,1),B(i)end do do i=1,m A(i,2)=1.0 end do if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E!invert(D)*E else if(m=n)then X=A.ix.B else write(*,*)No roots!end if open(unit=20,file=parameters.out)do i=1,n write(20,*)X(i)end do open(unit=30,file=data.out)write(30,(4a12)x,yexp,ycal,yserr%yerr_ave=0.yerr_max=0.do i=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)then yerr_max=(ycal-B(i)/B(i)*100 end if write(30,(f12.2,2f12.4,f12.3)A(i,1),B(i),ycal,(ycal-B(i)/B(i)*100 end do yerr_ave=yerr_ave/m write(30,*)write(30,*)write(30,(a20,f8.3)yerr_ave%=,yerr_ave write(30,(a20,f8.3)yerr_max%=,yerr_max write(30,*)close(30)stopend program least_squareRline函数函数n n拟合直线y=ax+b直线n nSubroutine RLINE(NOBS,XDATA,YDATA,B0,B1,STAT)integer NOBS:数据点数目 real XDATA(NOBS)：数据点X轴值 real yDATA(NOBS)：数据点Y轴值 real B0,B1:Y=B0+B1*X real stat(12):返回的信息。例：例：program main use IMSL implicit none integer,parameter:NOBS=5 real XDATA(NOBS),YDATA(NOBS)real B0,B1 real STAT(12)real xinc,xp,value real F,X F(X)=2*X+3 integer I!产生数据点产生数据点 do I=1,NOBS XDATA(I)=real(I)YDATA(I)=F(XDATA(I)end do call RLINE(NOBS,XDATA,YDATA,B0,B1,STAT)!结果一定为结果一定为y=2X+3,因为数据点是根据这个函数来产生的因为数据点是根据这个函数来产生的.write(*,(Y=,F5.2,X+F5.2)B1,B0 stopend programRcurv函数函数n n用最小二乘法来求一条最接近数据点的n次多项式。n nSubroutine RCURV(NOBS,XDATA,YDATA,NDEG,B,SSPOLY,STAT)integer NOBS:数据点数目 real XDATA(NOBS)：数据点X轴值 real yDATA(NOBS)：数据点Y轴值 integer NDEG:所要计算的多项式次数n nReal B(NDEG+1):多项式系数 real SSPOLY:返回信息 real STAT(10):返回信息例：例：program main use IMSL implicit none integer,parameter:NOBS=10,NDEG=2 real XDATA(NOBS),YDATA(NOBS)real B(NDEG+1),SSPOLY(NDEG+1)real STAT(12)real F,X,R F(X,R)=R+1+2*X+3*X*X integer Icall random_seed()!产生数据点产生数据点 do I=1,NOBS XDATA(I)=real(I)call random_number(R)YDATA(I)=F(XDATA(I),R)end do call RCURV(NOBS,XDATA,YDATA,NDEG,B,SSPOLY,STAT)write(*,(F5.2+F5.2X+F5.2X*X)B stopend program小小 结结n nF(x)=0求根 二分法，牛顿迭代法线性方程(n*n)组求解：高斯消元法，逆矩阵法超定线性方程(m*n,mn)组求解：最小二乘法性线和非线性拟合：最小二乘法应用举例应用举例n n如何求范德华方程中如何求范德华方程中气相气相的体积？的体积？n n方程展开以后为：方程展开以后为：对于水来说：对于水来说：应用程序应用程序program mainimplicit none real :T,P,V!units:K,bar,cm3/mol real :a,b real :V1,V2,V3 real :Ps real,parameter:R=83.14472!unit:cm3.bar/(mol.K)real,parameter:Tc=647.096 !unit:K real,parameter:Pc=220.64 !unit:bar100 write(*,*)write(*,*)Input T(K)and P(bar):read(*,*)T,P if(T=273.16.and.T=273.16.and.TPs)then write(*,*)Not vapor phase!goto 100 end if a=27.0*R*2*Tc*2/(64.0*Pc)b=R*Tc/(8.0*Pc)call roots(-(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(*,*)stopend program main!from Wagner and Pruss,J.Phys.Chem.Ref.Data,22(3),1993.subroutine H2Osaturation(Ts,Ps)!units:K,bar implicit none real :Ts,t,Ps real :a1,a2,a3,a4,a5,a6 real,parameter:Tc=647.096 !unit:K real,parameter:Pc=220.64 !unit:bar t=1.0-Ts/Tc a1=-7.85951783 a2=1.84408259 a3=-11.7866497 a4=22.6807411 a5=-15.9618719 a6=1.80122502 ps=Pc*exp(Tc/Ts*(a1*t+a2*t*1.5+a3*t*3.0+&a4*t*3.5+a5*t*4.0+a6*t*7.5)returnend subroutine H2Osaturationsubroutine roots(a,b,c,Z1,Z2,Z3)implicit none real :a,b,c real :p,q,r real :delta real :angle real :x1,x2,x3 !xi is the root of x3+px+q=0;real :z1,z2,z3,zmax !zi is the root of z3+az2+bz+c=0 real,external :power real,parameter:pi=3.14159265p=b-a*2/3.0q=2.0/27.0*a*3+c-a*b/3.0 delta=q*2/4.0+p*3/27.0if(delta 0)then x1=power(-q/2.0+sqrt(delta),1.0/3.0)+power(-q/2.0-sqrt(delta),1.0/3.0)x2=0 x3=0 Z1=x1-a/3.0 Z2=0 Z3=0 end if if(delta=0)then x1=2.0*power(-q/2.0,1.0/3.0)x2=power(q/2.0,1.0/3.0)x3=x2 Z1=x1-a/3.0 Z2=x2-a/3.0 Z3=Z2end ifif(delta=0)then power=x*y else power=-(-x)*y end if return end function power思思 考考n n如果采用牛顿迭代法，应用程序如何编写？如果采用牛顿迭代法，应用程序如何编写？n n这时初始体积可以取理想状态方程体积这时初始体积可以取理想状态方程体积n n即：即：