第一章客观分析

第1章客观分析1.1拉格朗日插值1.1.1FORTRAN程序LAGINT.F!!根据给定点的函数值，应用拉格朗日插值公式，计算指定插值点处的！函数值。！在指定点的前后各取四个结点，用七次拉格朗日插值公式。！X，Y为输入参数，各有N个点！T为指定插入点！Z为输出参数，是插值点T处的近似值。SUBROUTIONELAGINT（X，Y，N，T，Z）DIMENSIONX（N），Y（N）DOUBLEPRECISIONX，Y，T，X，SZ=0.0IF(N.LE.0)RETURNIF(N.EQ.1)THENZ=Y(1)ENDIFIF(N.EQ.2)THENZ=(Y(1)*(T-X(2))-Y(2)*(T-X(1)))/(X(1)-X(2))RETURNENDIFI=110IF(X(1).LT.T)I=I+1IF(I.LE.N)GOTO10ENDIFK=I-4２IF(K.LT.1)K=1M=I+3IF(M.GT.N)M=NDO30I=K,MS=1.0DO20J=K,MIF(J.NE.I)THENS=S*(T-X(J))/(X(I)-X(J))ENDIF20CONTINUEZ=Z+S*Y(I)30CONTINUERETURNEND1.2样条插值FORTRAN程序SPLINE.F!SPLINE.FOR--NaturalCubicSpline!CallSplineCalculatewiththetwoarraysofcontrolpoints(and!thenumberofelementsineacharray).Tointerpolatethespline,!CallSplineEvaluateforeach'X'forwhichyouwanta'Y'alongthe!spline'scurve.!subroutineSplineCalculate(inx,ina,numx)!ms$if.not.defined(LINKDIRECT)!ms$attributesdllexport::SplineCalculate!ms$endifimplicitnone３real*4ina(*),inx(*)integernumxintegerMAXPOINTSparameter(MAXPOINTS=1000)integernxsreal*4x(0:MAXPOINTS)real*4a(0:MAXPOINTS),b(0:MAXPOINTS)real*4c(0:MAXPOINTS),d(0:MAXPOINTS)common/SplineData/a,b,c,d,nxs,xintegerireal*4alpha(0:MAXPOINTS),l(0:MAXPOINTS)real*4mu(0:MAXPOINTS),z(0:MAXPOINTS)real*4h(0:MAXPOINTS)!=====Step1nxs=numx-1doi=0,nxsx(i)=inx(i+1)a(i)=ina(i+1)enddodoi=0,nxs-1h(i)=x(i+1)-x(i)enddox(nxs+1)=40000!=====Step2doi=1,nxs-1alpha(i)=3.0*(a(i+1)*h(i-1)-a(i)*４&(x(i+1)-x(i-1))+a(i-1)*h(i))/&(h(i-1)*h(i))enddo!=====Step3l(0)=1.0mu(0)=0.0z(0)=0.0!=====Step4doi=1,nxs-1l(i)=2.0*(x(i+1)-x(i-1))-h(i-1)*mu(i-1)mu(i)=h(i)/l(i)z(i)=(alpha(i)-h(i-1)*z(i-1))/l(i)enddo!=====Step5l(nxs)=1.0z(nxs)=0.0c(nxs)=0.0!=====Step6doi=nxs-1,0,-1c(i)=z(i)-mu(i)*c(i+1)b(i)=(a(i+1)-a(i))/h(i)-h(i)*&(c(i+1)+2.0*c(i))/3.0d(i)=(c(i+1)-c(i))/(3.0*h(i))enddoend５subroutineSplineEvaluate(evalx,evaly)!ms$if.not.defined(LINKDIRECT)!ms$attributesdllexport::SplineEvaluate!ms$endifreal*4evalx,evalyintegerMAXPOINTSparameter(MAXPOINTS=1000)integernxsreal*4x(0:MAXPOINTS)real*4a(0:MAXPOINTS),b(0:MAXPOINTS)real*4c(0:MAXPOINTS),d(0:MAXPOINTS)common/SplineData/a,b,c,d,nxs,xreal*4termintegerii=0dowhile(.not.((x(i).le.evalx).and.(x(i+1).gt.evalx)))i=i+1enddoterm=evalx-x(i)evaly=a(i)+b(i)*term+c(i)*term*term+&d(i)*term*term*termend６1.3曲线拟合1.3.1.1FORTRAN程序!HPINT.FOR!用最小二乘法求N个数据点的拟合多项式。！X，Y输入参数，分别存放N个数据点的X坐标与Y坐标！A输出参数，存放M-1次拟合多项式的M个系数！N输入参数，数据点数！M输入参数，拟合多项式的项数，MN且M20！DT1，DT2，DT3—输出参数，分别为拟合多项式与数据点偏差的平方和、绝对！值之和、绝对值的最大值！SUBROTINEHPINT(X,Y,A,N,M,DT1,DT2,DT3)DIMENSIONX(N),Y(N),A(M),S(20),T(20),B(20)DOUBLEPRECISIONX,Y,A,S,T,B,DT1,DT2,DT3,Z,D1,P,&C,D2,G,Q,DTDO5I=1,M5A(I)=0.0IF(M.GT.N)M=NIF(M.GT.20)M=20Z=0.0DO10I=1,N10Z=Z+X(I)/NB(1)=1.0D1=NP=0.0C=0.0DO20I=1,NP=P+(X(I)–Z)C=C+Y(I)７20CONTINUEC=C/D1P=P/D1A(1)=C*B(1)IF(M.GT.1)THENT(2)=1.0T(1)=-PD2=0.0C=0.0G=0.0DO30I=1,NQ=X(I)–Z–PD2=D2+Q*QC=Y(I)*Q+CG=(X(I)–Z)*Q*Q+G30CONTINUEC=C/D2P=G/D2Q=D2/D1D1=D2A(2)=C*T(2)A(1)=C*T(1)+A(1)ENDIFDO100J=3,MS(J)=T(J-1)S(J-1)=-P*T(J-1)+T(J-2)IF(J.GE.4)THENDO40K=J-2,2,-140S(K)=-P*T(K)+T(K-1)–Q*B(K)ENDIF８S(1)=–P*T(1)–Q*B(1)D2=0.0C=0.0G=0.0DO70I=1,NQ=S(J)DO60K=J-1,1,-160Q=Q*(x(I)-z_+s(k)D2=D2+Q*QC=Y(I)*Q+CG=(X(I)-Z)QQ+G70CONTINUECC/D2P=G/D2Q=D2/D1D1=D2A(J)=C*S(J)T(J)=S(J)DO80K=J-1,1,-1A(K)=C*S(K)+A(K)B(K)=T(K)T(K)=S(K)80CONTINUE100CONTINUEDT1=0.0DT2=0.0DT3=0.0DO120I=1,NQ=A(M)DO110K=M-1,1,-1９110Q=Q*(X(I)-Z)+A(K)DT=Q-Y(I)IF(ABS(DT).GT.DT3)D3=ABS(DT)DT1=DT1+DT*DTDT2=DT2+ABS(DT)120CONTINUERETURNEND