!****************************************************************** !* Program to demonstrate Lagrange derivative interpolation * !* -------------------------------------------------------------- * !* Ref.: Basic Scientific Subroutines Vol. II page 314 * !* By F.R. Ruckdeschel. BYTE/McGRAW-HILL, 1981 [BIBLI 01] * !* * !* F90 version by J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* -------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* X 2COS(2X) YP YP1 ERROR 1 ERROR 2 * !* -------------------------------------------------------------- * !* 0.00 2.000000 1.999961 2.006643 -0.0000394 0.0066434 * !* 0.05 1.990008 1.989970 1.996569 -0.0000385 0.0065603 * !* 0.10 1.960133 1.960096 1.966545 -0.0000373 0.0064118 * !* 0.15 1.910673 1.910637 1.916872 -0.0000357 0.0061991 * !* 0.20 1.842122 1.842088 1.848047 -0.0000337 0.0059246 * !* 0.25 1.755165 1.755134 1.760756 -0.0000314 0.0055908 * !* 0.30 1.650671 1.650642 1.655872 -0.0000288 0.0052011 * !* 0.35 1.529684 1.529659 1.534444 -0.0000259 0.0047595 * !* 0.40 1.393413 1.393391 1.397684 -0.0000227 0.0042704 * !* 0.45 1.243220 1.243201 1.246958 -0.0000193 0.0037386 * !* 0.50 1.080605 1.080589 1.083774 -0.0000157 0.0031694 * !* 0.55 0.907192 0.907180 0.909761 -0.0000120 0.0025685 * !* 0.60 0.724715 0.724707 0.726658 -0.0000081 0.0019420 * !* 0.65 0.534998 0.534993 0.536294 -0.0000042 0.0012961 * !* 0.70 0.339934 0.339934 0.340571 -0.0000002 0.0006372 * !* 0.75 0.141474 0.141478 0.141446 0.0000038 -0.0000280 * !* 0.80 -0.058399 -0.058391 -0.059092 0.0000078 -0.0006929 * !* 0.85 -0.257689 -0.257677 -0.259040 0.0000116 -0.0013510 * !* 0.90 -0.454404 -0.454389 -0.456400 0.0000154 -0.0019955 * !* 0.95 -0.646579 -0.646560 -0.649199 0.0000190 -0.0026201 * !* 1.00 -0.832294 -0.832271 -0.835512 0.0000224 -0.0032185 * !* -------------------------------------------------------------- * !****************************************************************** PROGRAM Derivati parameter(NMAX=100,pas=0.05) real*8 X(NMAX),Y(NMAX) real*8 xx,yy,yy1 integer error,i,n,ndata !derivative of 2*sin(x)*cos(x) between 0 and 1 n=4 !level of Lagrange interpolation ndata=26 !number of table points ! building X & Y Tables do i = 1, ndata X(i) = pas*(i-1) Y(i) = 2.d0*cos(X(i))*sin(X(i)) end do xx=0.d0 print *,' X 2COS(2X) YP YP1 ERROR 1 ERROR 2 ' !heading print *,'--------------------------------------------------------------' !main loop of derivation do while (xx < 1.d0+pas) call Deriv(ndata,n,X,Y,xx,yy,error) call Deriv1(ndata,X,Y,xx,yy1,error) if (error.eq.0) then write(*,50) xx,2.d0*dcos(2.d0*xx),yy,yy1, & yy-2.d0*dcos(2.d0*xx), yy1-2.d0*dcos(2.d0*xx) end if xx=xx+pas end do print *,'--------------------------------------------------------------' 50 format(' ',f4.2,' ',f9.6,' ',f9.6,' ',f9.6,' ',f10.7,' ',f10.7) stop end Subroutine Deriv(n, nl, X, Y, x1, Yp, error) !***************************************************** !* Lagrange derivative interpolation procedure Deriv * !* NL is the level of the interpolation ( ex. NL=3 ) * !* N is the total number of table values. * !* X[i], Y[i] are the coordinate table values, Y[i] * !* being the dependant variable. The X[i] may be * !* arbitrarily spaced. x1 is the interpolation point * !* which is assumed to be in the interval with at * !* least one table value to the left, and NL to the * !* right. Yp is returned as the desired derivative. * !* error is set at TRUE if x1 is not in the interval.* !***************************************************** parameter(NMAX=100) integer n,nl,error real*8 X(NMAX),Y(NMAX) real*8 x1,Yp integer i,j,k,ll real*8 L(10), M(10,10) error=0 ! x1 not in interval [1:N-NL] if (x1X(n-nl)) then error=1 print *,' STOP: x not between X[1] or X[N-4].' end if if (error.eq.0) then i=0; do while (x1 >= X(i)) i=i+1 end do i=i-1 do j = 0, nl L(j)=0.d0 do k = 0, nl M(j,k)=1.d0 end do end do Yp=0.d0 do k = 0, nl do j = 0, nl if (j.ne.k) then do ll = 0, nl if (ll.ne.k) then if (ll.eq.j) then M(ll,k)=M(ll,k)/(X(i+k)-X(i+j)) else M(ll,k)=M(ll,k)*(x1-X(j+i))/(X(i+k)-X(i+j)) end if end if end do end if end do do ll = 0, nl if (ll.ne.k) L(k)=L(k)+M(ll,k) end do Yp=Yp+L(k)*Y(i+k) end do end if return end ! Deriv() !********************************************************* !CALCULE LES COEFFICIENTS A,B DE LA PARABOLE Y=A*X*X+B*X+C !PASSANT PAR LES 3 POINTS : (X1,Y1),(X2,Y2) ET (X3,Y3) !COEFFICIENT C NON UTILISE ICI. !--------------------------------------------------------- !Calculates coefficients, a, b of parabola Y=A*X+X+B*X=C !passing through 3 points: (X1,Y1), (X2,Y2) and (X3,Y3). !Coefficient c is not used here. !********************************************************* Subroutine PARABOLE(x1,y1,x2,y2,x3,y3,a,b) real*8 x1,y1,x2,y2,x3,y3,a,b real*8 alpha,beta,gamma,delta alpha=x1*x1-x2*x2 beta=x1-x2 gamma=x2*x2-x3*x3 delta=x2-x3 a=(y1-2.d0*y2+y3)/(alpha-gamma) b=(y1-y2-alpha*a)/beta return end !**************************************************** !* Interpolation of order=2 ( parabola ) * !* by J-P Moreau * !**************************************************** Subroutine Deriv1(n, X, Y, x1, Yp, error) parameter(NMAX=100) real*8 X(NMAX),Y(NMAX) real*8 a,b,x1,Yp integer error,i,n error=0 if (x1X(n-2)) then error=1 return end if i=0 do while (x1>=X(i)) i=i+1 end do i=i-1 call PARABOLE(X(i),Y(i),X(i+1),Y(i+1),X(i+2),Y(i+2),a,b) !Derivative in x1 Yp=2.d0*a*x1 + b return end ! End of file derivati.f90