MODULE DEFINE_FCN ! Declare global parameters. DOUBLE PRECISION, PARAMETER :: Q=5D0,& PI=3.141592653589793D0 ! Define global variables for the number of differential ! equations, NODE, the number of unknown parameters, NPAR, ! and the number of boundary conditions at the left end of ! the interval, LEFTBC. The number of boundary conditions ! at the right end is NODE+NPAR-LEFTBC. INTEGER, PARAMETER :: NODE=2,NPAR=1,& LEFTBC=2,RIGHTBC=NODE+NPAR-LEFTBC CONTAINS SUBROUTINE FSUB(X,Y,P,DYDX) DOUBLE PRECISION :: X,Y(NODE),P(NPAR),LAMBDA,DYDX(NODE) LAMBDA = P(1) DYDX(1) = Y(2) DYDX(2) = -(LAMBDA - 2D0*Q*COS(2D0*X))*Y(1) END SUBROUTINE FSUB SUBROUTINE DFSUB(X,Y,P,DFDY,DFDP) DOUBLE PRECISION :: X,Y(NODE),P(NPAR),LAMBDA,& DFDY(NODE,NODE),DFDP(NODE,NPAR) LAMBDA = P(1) DFDY = 0D0 DFDY(1,2) = 1D0 DFDY(2,1) = -(LAMBDA - 2D0*Q*COS(2D0*X)) DFDP = 0D0 DFDP(2,1) = - Y(1) END SUBROUTINE DFSUB SUBROUTINE BCSUB(YA,YB,P,BCA,BCB) DOUBLE PRECISION :: YA(NODE),YB(NODE),P(NPAR),& BCA(LEFTBC),BCB(RIGHTBC) BCA(1) = YA(2) BCA(2) = YA(1) - 1D0 BCB(1) = YB(2) END SUBROUTINE BCSUB SUBROUTINE DGSUB(YA,YB,DYA,DYB,P,DAP,DBP) DOUBLE PRECISION :: YA(NODE),YB(NODE),& DYA(LEFTBC,NODE),DYB(RIGHTBC,NODE),& P(NPAR),DAP(LEFTBC,NPAR),DBP(RIGHTBC,NPAR) DYA = 0D0 DYA(1,2) = 1D0 DYA(2,1) = 1D0 DYB = 0D0 DYB(1,2) = 1D0 DAP = 0D0 DBP = 0D0 END SUBROUTINE DGSUB SUBROUTINE GUESS_Y(X,Y) DOUBLE PRECISION :: X,Y(NODE) Y(1) = COS(4D0*X) Y(2) = -4D0*SIN(4D0*X) END SUBROUTINE GUESS_Y END MODULE DEFINE_FCN !************************************************************ PROGRAM MATHIEU ! This is the MAT4BVP example of Matlab. It finds the fourth ! eigenvalue of Mathieu's equation ! ! y'' + (lambda - 2*q*cos(2*x))*y = 0 ! ! on the interval [0, pi] with boundary conditions y'(0) = 0, ! y'(pi) = 0 when the parameter q = 5. ! ! Special codes for Sturm-Liouville problems can compute specific ! eigenvalues. The general-purpose code BVP_SOLVER can only compute ! an eigenvalue near to a guessed value. We can make it much more ! likely that we compute the eigenfunction corresponding to the fourth ! eigenvalue by supplying a guess that has the correct qualitative ! behavior. The eigenfunction y(x) is determined only to a constant ! multiple, so the normalizing condition y(0) = 1 is used to specify ! a particular solution. USE DEFINE_FCN USE BVP_M TYPE(BVP_SOL) :: SOL ! Working variables for output. INTEGER :: I DOUBLE PRECISION :: XPLOT,YPLOT(NODE),LAMBDA(NPAR) DOUBLE PRECISION :: A=0D0,B=PI ! BVP_INIT is used to form an initial guess for a mesh of 10 equally spaced ! points. The guess cos(4x) for y(x) and its derivative as guess for y'(x) ! are evaluated in GUESS_Y. The desired eigenvalue is the one nearest the ! guess lambda = 15. A guess for unknown parameters is supplied in an array ! as the last argument of BVP_INIT. SOL = BVP_INIT(NODE,LEFTBC,(/ A,B /),GUESS_Y,(/ 15D0 /)) SOL = BVP_SOLVER(SOL,FSUB,BCSUB,DFDY=DFSUB,DBCDY=DGSUB) ! Call BVP_EVAL to get the approximation for the parameter. ! Parameters are held in an array even when there is only ! one unknown parameter, LAMBDA(1). The parameters are also ! available as SOL%PARAMETERS. CALL BVP_EVAL(SOL,LAMBDA) PRINT *,'Computed the parameter lambda to be ', LAMBDA(1) PRINT *,'The mat4bvp program of MATLAB says it is about 17.097.' PRINT *,' ' ! Evaluate the solution and write it to an output file. OPEN(UNIT=7,FILE="mathieu.dat") WRITE(UNIT=7,FMT="(2D14.5)") LAMBDA(1),LAMBDA(1) DO I = 1,100 XPLOT = A + (I-1)*(B - A)/99D0 CALL BVP_EVAL(SOL,XPLOT,YPLOT) WRITE(UNIT=7,FMT="(2D12.4)") XPLOT,YPLOT(1) END DO PRINT *,'The solution evaluated at 100 equally spaced points' PRINT *,'can be imported into Matlab with' PRINT *,' ' PRINT *,' >> [x,y] = mathieu;' PRINT *,' ' PRINT *,'For this the output file "mathieu.dat" must be' PRINT *,'in the same directory as mathieu.m.' PRINT *,' ' PRINT *,'After importing the data, mathieu.m plots the solution.' PRINT *,' ' CALL BVP_TERMINATE(SOL) END PROGRAM MATHIEU