最小二乘法应用幻灯片.ppt

《最小二乘法应用幻灯片.ppt》由会员分享，可在线阅读，更多相关《最小二乘法应用幻灯片.ppt（85页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、最小二乘法应用最小二乘法应用第1页，共85页，编辑于2022年，星期六主要内容主要内容n n1.基本原理n n2.程序编写n n3.应用举例第2页，共85页，编辑于2022年，星期六AXAX=b b一般一般m*n线性方程组线性方程组(mn)1.基本原理基本原理第3页，共85页，编辑于2022年，星期六超定方程组超定方程组：方程的个数超过未知量的个数，一般来说：方程的个数超过未知量的个数，一般来说是没有解的。即对于任意一组数是没有解的。即对于任意一组数(x1,x2,xn)不会全为零。我们设想去求一组数不会全为零。我们设想去求一组数第4页，共85页，编辑于2022年，星期六使得：使得：取得最小值取

2、得最小值利用多元函数求极值的方法可得下式利用多元函数求极值的方法可得下式第5页，共85页，编辑于2022年，星期六即：即：用矩阵形式给出，即第6页，共85页，编辑于2022年，星期六n n可以证明如果如果A是列满秩的，则该方程存在唯一解，且使得取得最小值，该解就称为超定方程组的最小二乘解。第7页，共85页，编辑于2022年，星期六例例 1：用最小二乘法求下列超定用最小二乘法求下列超定方程组的近似解方程组的近似解第8页，共85页，编辑于2022年，星期六n n解：第9页，共85页，编辑于2022年，星期六故得方程组故得方程组解得解得 x1=0.79272,x2=-1.4641第10页，共85页，

3、编辑于2022年，星期六2.程序编写程序编写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/第11页，共85页，编辑于2022年，星期六 if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X

4、=D.ix.E!invert(D)*E else if(m=n)then X=A.ix.B else write(*,*)No roots!end if 第12页，共85页，编辑于2022年，星期六 do i=1,n write(*,*)X(i)end do stopend program least_square第13页，共85页，编辑于2022年，星期六如何拟合线性如何拟合线性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？第14页，共85页，编辑于2022年，星期六prog

5、ram 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/第15页，共85页，编辑于2022年，星期六if(m n)then C=transpose(a)D=matmul(c,a)E=matmul(c,b)X=D.ix.E

6、!invert(D)*Eelse if(m=n)then X=A.ix.B else write(*,*)No roots!end if第16页，共85页，编辑于2022年，星期六 do i=1,n write(*,*)X(i)end do stopend program least_square第17页，共85页，编辑于2022年，星期六n nPs.in3.应用举例应用举例:O2饱和压的拟合与计算饱和压的拟合与计算550.001785755.50.0020802560.00241656.50.002798570.003231457.50.0037217580.00427558.50.0048

7、98590.005597759.50.0063817600.007258260.50.008236610.009324361.50.010533620.01187362.50.013356630.01499363.50.016797640.01878164.50.020961第一列为温度，单位为第一列为温度，单位为K；第二列为实验的饱和压，单位为第二列为实验的饱和压，单位为bar.第18页，共85页，编辑于2022年，星期六650.02334965.50.025964660.02881966.50.031934670.03532667.50.039015680.04301968.50.0473

8、6690.05205969.50.057139700.06262370.50.068535710.07490171.50.081746720.08909872.50.096984730.1054373.50.11448740.1241474.50.13446第19页，共85页，编辑于2022年，星期六750.1454775.50.15721760.1696976.50.18297770.1970877.50.21205780.2279278.50.24474790.2625379.50.28135800.3012380.50.32222810.3443681.50.36769820.39226

9、82.50.41812830.445383.50.47386840.5038584.50.53532第20页，共85页，编辑于2022年，星期六850.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.5559第21页，共85页，编辑于2022年，星期六951.630895.51.708

10、5961.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.6438第22页，共85页，编辑于2022年，星期六1053.7853105.53.93091064.0806106.54.23441074.3925107.54.55491084.7217108.54.89291095.0687109.55.2491105.4341

11、10.55.62371115.8183111.56.01761126.222112.56.43131136.6458113.56.86541147.0902114.57.3204第23页，共85页，编辑于2022年，星期六1157.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.7981

12、24.513.15第24页，共85页，编辑于2022年，星期六12513.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.737第25页，共85页，编辑于2022年，星期六13522.25135.522.77313623.303136.523.84313724.392

13、137.524.9513825.516138.526.09313926.678139.527.27314027.878140.528.49214129.116141.529.7514230.394142.531.04914331.713143.532.38814433.074144.533.77第26页，共85页，编辑于2022年，星期六14534.477145.535.19614635.925146.536.66614737.418147.538.18214838.958148.539.74614940.547149.541.3615042.186150.543.02515143.87815

14、1.544.74515245.626152.546.52215347.434153.548.36215449.307154.550.271第27页，共85页，编辑于2022年，星期六要求：要求：n n实验数据共200个，要求从文件（文件名：Ps.in)读入，拟合参数，计算出饱和压与实验数据间的误差，找出最大误差和平均误差，并以文件（文件名：Ps.out）形式输出。n n温度、饱和压及相应中间变量都采取双精度实数。第28页，共85页，编辑于2022年，星期六氧气饱合压拟合公式如下：氧气饱合压拟合公式如下：其中其中：c(1)=？c(2)=？编程计算c(3)=？c(4)=？Tc=154.581!un

15、it:k Pc=50.43 !unit:bar第29页，共85页，编辑于2022年，星期六输出文件输出文件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.

16、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

17、.50 0.0210 0.0210 0.054第30页，共85页，编辑于2022年，星期六 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

18、 49.3043 -0.005 154.50 50.2710 50.2708 0.000 Pserr_ave%=0.015 Pserr_max%=-0.157第31页，共85页，编辑于2022年，星期六程程 序序!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)

19、,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 第32页，共85页，编辑于2022年，星期六 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)we

20、ight(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)第33页，共85页，编辑于2022年，星期六 open(unit=30,file=Ps.out)write(30,(4a12)T(K),Psexp(bar),Psc

21、al,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 e

22、nd if write(30,(f12.2,2f12.4,f12.3)T(i),Ps(i),Pscal(i),(Pscal(i)-Ps(i)/Ps(i)*100 end do第34页，共85页，编辑于2022年，星期六 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)第35页，共85页，编辑于2022年，星期六!-to calculate Ps at a

23、ny 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 main第36页，共85页，编辑于2022年，星期六subr

24、outine 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第37页，共85页，编辑于2022年，星期六 !COE=0.0

25、!COE1=0.0 Do k=1,4 Do j=1,4 COE(k,j)=0.0 end do COE1(k)=0.0 End do第38页，共85页，编辑于2022年，星期六 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 xis

26、hu第39页，共85页，编辑于2022年，星期六!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 ENDIF第40页，共85页，编辑于2022年，星期六210 CONTINUE IF(D=1

27、d-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)=T第41页，共85页，编辑于2022年，星期六230 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,

28、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)第42页，共85页，编辑于2022年，星期六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第50

29、页，共85页，编辑于2022年，星期六求解超定非线性方程组？求解超定非线性方程组？n n要使用非线性拟合程序，比较复杂。通常有几个未知量，就要给几个初值。往往每给定一组初值，由目标函数得到的最小平均误差也不一样。所以要根据实际情况，进行多次尝试和分析，得出最后的参数值。第51页，共85页，编辑于2022年，星期六可化为多元线性回归可化为多元线性回归函数式 通式：变换变量变换常数yxabcylnxcosxcabxabcx/yxabcy1/xabclnyxlnabc第52页，共85页，编辑于2022年，星期六作业程序作业程序1!This program is to simulate paramet

30、ers 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第53页，共85页，编辑于2022年，星期六 open(unit=10,file=data.in)do i

31、=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)第54页，共85页，编辑于2022年，星期六 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

32、(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第55页，共85页，编辑于2022年，星期六 yerr_ave=yerr_ave/n write(30,*)write(30,*)write(30,(a20,f8.3)yerr_ave%=,yerr_ave write(30,

33、(a20,f8.3)yerr_max%=,yerr_max write(30,*)close(30)stopend program main第56页，共85页，编辑于2022年，星期六subroutine 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)c

34、oe1=matmul(b,y)returnend subroutine xishu第57页，共85页，编辑于2022年，星期六作业程序作业程序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.

35、3,19.1,20.8/第58页，共85页，编辑于2022年，星期六 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第59页，共85页，编辑于2022年，星期六do i=1,n write(*,*)X(i)end do stopend program least_square第60页，共85页，编辑于2022年，星期六作业程序作业程序3program least_square use

36、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第61页，共85页，编辑于2022年，星期六 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第62页，共85页，编辑于2022年

37、，星期六 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第63页，共85页，编辑于2022年，星期六 open(unit=30,file=data.out)write(30,(4a12)x,yexp,ycal,yserr%yerr_ave=0.yerr_max

38、=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第64页，共85页，编辑于2022年，星期六 yerr_ave=yerr_ave/m write(30,*)write(30,*)write(

39、30,(a20,f8.3)yerr_ave%=,yerr_ave write(30,(a20,f8.3)yerr_max%=,yerr_max write(30,*)close(30)stopend program least_square第65页，共85页，编辑于2022年，星期六Rline函数函数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*

40、X real stat(12):返回的信息。第66页，共85页，编辑于2022年，星期六例：例：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第67页，共85页，编辑于2022年，星期六!产生数据点产生数据点 do I=1,NOBS XDATA(I)=real(I)YDATA(I)=F(XDATA(I)end do cal

41、l RLINE(NOBS,XDATA,YDATA,B0,B1,STAT)!结果一定为结果一定为y=2X+3,因为数据点是根据这个函数来产生的因为数据点是根据这个函数来产生的.write(*,(Y=,F5.2,X+F5.2)B1,B0 stopend program第68页，共85页，编辑于2022年，星期六Rcurv函数函数n n用最小二乘法来求一条最接近数据点的n次多项式。n nSubroutine RCURV(NOBS,XDATA,YDATA,NDEG,B,SSPOLY,STAT)integer NOBS:数据点数目 real XDATA(NOBS)：数据点X轴值 real yDATA(N

42、OBS)：数据点Y轴值 integer NDEG:所要计算的多项式次数第69页，共85页，编辑于2022年，星期六n nReal B(NDEG+1):多项式系数 real SSPOLY:返回信息 real STAT(10):返回信息第70页，共85页，编辑于2022年，星期六例：例：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

43、+1+2*X+3*X*X integer I第71页，共85页，编辑于2022年，星期六call 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第72页，共85页，编辑于2022年，星期六小小 结结n nF(x)=0求根 二分法，牛顿迭代法线性方程(n*n

44、)组求解：高斯消元法，逆矩阵法超定线性方程(m*n,mn)组求解：最小二乘法性线和非线性拟合：最小二乘法第73页，共85页，编辑于2022年，星期六应用举例应用举例n n如何求范德华方程中气相的体积？第74页，共85页，编辑于2022年，星期六n n方程展开以后为：对于水来说：对于水来说：第75页，共85页，编辑于2022年，星期六应用程序应用程序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:cm

45、3.bar/(mol.K)real,parameter:Tc=647.096 !unit:K real,parameter:Pc=220.64 !unit:bar第76页，共85页，编辑于2022年，星期六100 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第77页，共85页，编辑于2022年，星期六 a=27.0*R*2*Tc*2/(64.0*Pc)b=R*Tc/(8.0*

46、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第78页，共85页，编辑于2022年，星期六!from Wagner and Pruss,J.Phys.Chem.Ref.Data,22(3),1993.subroutine H2Osaturation(Ts,Ps)!units:K,bar implicit none real :Ts,t

47、,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 H2Osaturation第7

48、9页，共85页，编辑于2022年，星期六subroutine 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.14159265第80页，共85页，编辑于2022年，星期六p=b-a*2/3.0q=2.0/2

49、7.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 第81页，共85页，编辑于2022年，星期六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 if第82页，共85页，编辑于2022年，星期六if(delta=0)then power=x*y else power=-(-x)*y end if return end function power第84页，共85页，编辑于2022年，星期六思思 考考n n如果采用牛顿迭代法，应用程序如何编写？n n这时初始体积可以取理想状态方程体积n n即：第85页，共85页，编辑于2022年，星期六