4: DOCUMENTATION MUSN **************** SPECIFICATION MUSN **************** SUBROUTINE MUSN(FDIF,Y0T,G,N,A,B,ER,TI,NTI,NRTI,AMP,ITLIM,Y,Q,U, 1 NU,D,PHI,KP,W,LW,IW,LIW,WG,LWG,IERROR) C C DOUBLE PRECISION A,B,ER(5),TI(NTI),AMP,Y(N,NTI),Q(N,N,NTI), C 1 U(NU,NTI),D(N,NTI),PHI(NU,NTI),W(LW),WG(LWG) C INTEGER N,NTI,NRTI,ITLIM,NU,KP,LW,IW(LIW),LIW,LWG,IERROR C EXTERNAL FDIF,Y0T,G **************** Purpose **************** MUSN solves the nonlinear two-point BVP dy(t)/dt = f(t,y) A <= t <= B or B <= t <= A g(y(a),y(b)) = 0 , where y(t) and f(t,y) are N-vector functions. **************** Method **************** MUSN uses a multiple shooting method for computing an approximate solution of the BVP at specified output points, which are also used as shooting points. If necessary, more output points (shooting points) are inserted during computa- tion. For integration a fixed grid is used. Output points are also grid points and the minimum number of grid points between two output points is 5. **************** Parameters **************** FDIF SUBROUTINE, supplied by the user with specification: SUBROUTINE FDIF(T,Y,F) DOUBLE PRECISION T,Y(N),F(N) where N is the order of the system. FDIF must evaluate the righthand- side of the differential equation, f(t,y) for t=T and y=Y and places the result in F(1),...,F(N). FDIF must be declared as EXTERNAL in the (sub)program from which MUSN is called. Y0T SUBROUTINE, supplied by the user with specification SUBROUTINE Y0T(T,Y) DOUBLE PRECISION T,Y(N) where N is the order of the system. Y0T must evaluate the initial appro- ximation y0(t) of the solution, for any value t=T and place the result in Y(1),...,Y(N). Y0T must be declared as EXTERNAL in the (sub)program from which MUSN is called. G SUBROUTINE, supplied by the user with specification SUBROUTINE G(N,YA,YB,FG,DGA,DGB) DOUBLE PRECISION YA(N),YB(N),FG(N),DGA(N,N),DGB(N,N) where N is the order of the system. G must evaluate g(y(A),y(B)) for y(A)=YA and y(B)=YB and place the result in FG(1),...,FG(N). Moreover G must evaluate the Jacobians dg(u,v)/du for u = YA and dg(u,v)/dv for v = YB and place the result in the arrays DGA anb DGB respectively. G must be declared as EXTERNAL in the (sub)program from which MUSN is called. N INTEGER, the order of the system. Unchanged on exit. A,B DOUBLE PRECISION, the two boundary points. Unchanged on exit. ER DOUBLE PRECISION array of dimension (5). On entry: ER(1) must contain the required tolerance for solving the differential equation. ER(2) must contain the initial tolerance with which a first aproximate solution will be computed. This approximation is then used as an initial approximation for the computation of an solution with an tolerance ER(2)*ER(2) and so on until the required tolerance is reached. As an initial tolerance max(ER(1),min(ER(2),1.d-2)) will be used. ER(3) must contain the machine precision. On exit: ER(1), ER(2) and ER(3) are unchanged. ER(4) contains an estimation of the condition number of the BVP. ER(5) contains an estimated error amplification factor. TI DOUBLE PRECISION array of dimension (NTI). On entry : if NRTI = 1 TI must contain the output points in strict monotone order: A=TI(1) < TI(2) <...< TI(n)=B. On exit TI(j), j=1,...,NRTI contains the output points. NTI INTEGER, NTI is one of the dimension of TI, X, S, Q, U en PHI. NTI must be greater than or equal to the total number of necessary output points + 1 (i.e. if the entry value for NRTI > 1, NTI may be equal to the entry value of NRTI + 1). Unchanged on exit. NRTI INTEGER. On entry NRTI is used to specify the output points. There are 3 ways to specify the output points: 1) NRTI = 0, the output points are determined automatically using AMP. 2) NRTI = 1, the output points are supplied by the user in the array TI. 3) NRTI > 1, the subroutine computes the (NRTI+1) output points TI(k) by TI(k) = A + (k-1) * (B - A) / NRTI ; so TI(1) = A and TI(NRTI+1) = B. Depending on the allowed increment between two succesive output points, more output points may be inserted in cases 2 and 3. On exit NRTI contains the total number of output points. AMP DOUBLE PRECISION. On entry AMP must contain the allowed increment between two output points. AMP is used to determine output points and to assure that the increment between two output points is at most AMP*AMP. A small value for AMP may result in a large number of output points. Unless 1 < AMP < .25 * sqrt(ER(1)/ER(3)) the default value .25 * sqrt(ER(1)/ER(3)) is used. Unchanged on exit. ITLIM INTEGER, maximum number of allowed iteration. Y DOUBLE PRECISION array of dimension (N,NTI). On exit Y(.,i), i=1,...,NRTI contains the solution at the output points TI(i), i=1,...,NRTI. Q DOUBLE PRECISION array of dimension (N,N,NTI). On exit Q(.,.,i), i=1,...,NRTI contains the orthogonal factors of the incremental recursion. U DOUBLE PRECISION array of dimension (NU,NTI). On exit U(.,i), i=2,...,NRTI contains the upper triangular factors of the incremental recursion. The elements are stored column wise, the j-th column of U is stored in U(nj+1,.),U(nj+2,.),...,U(nj+j,.), where nj = (j-1) * j / 2. NU INTEGER, NU is one of the dimension of U and PHI. NU must be greater than or equal to N * (N+1) / 2. Unchanged on exit. D DOUBLE PRECISION array of dimension (N,NTI). On exit D(.,i) i=2,...,NRTI contain the inhomogeneous terms of the incremental recursion. PHI DOUBLE PRECISION array of dimension (NU,NTI). On exit PHI(.,i), i=1,...,NRTI contains the fundamental solution of the incremental recursion. The fundamental solution is upper triangular and stored in the same way as the upper triangular U. KP INTEGER, On exit KP contains the dimension of the increasing solution space. W DOUBLE PRECISION array of dimension (LW). Used as work space. LW INTEGER. LW is the dimension of W. LW >= 7*N + 3*N*NTI + 4*N*N . IW INTEGER array of dimension (LIW). Used as work space. LIW INTEGER. LIW is the dimension of IW. LIW >= 3*N + NTI . WG DOUBLE PRECISION array of dimension (LWG). WG is used to store the integration grid points. LWG INTEGER. LWG is the dimension of WG. LWG >= (total number of grid points) / 5. The minimum number of grid points between 2 succesive output points is 5, so the minimum value for LWG is the number of actually used output points. Initially a crude estimate for LWG has to be made (see also IERROR 219). IERROR INTEGER, error indicator. On entry : if IERROR=1, diagnostics will be printed during computation. On exit : if IERROR = 0 no errors have been detected. **************** Error indicators **************** IERROR 0 No errors detected. 101 INPUT error: either ER(1) < 0 or ER(2) < 0 or ER(3) < 0. TERMINAL ERROR. 105 INPUT error: either N < 1 or NRTI < 0 or NTI < 3 or NU < N*(N+1)/2 or A=B. TERMINAL ERROR. 106 INPUT error: either LW < 7*N + 3*N*NTI + 4*N*N or LIW < 3*N + NTI. TERMINAL ERROR. 120 INPUT error: the routine was called with NRTI=1, but the given output points in the array TI are not in strict monotone order. TERMINAL ERROR. 121 INPUT error: the routine was called with NRTI=1, but the first given output point or the last output point is not equal to A or B. TERMINAL ERROR. 122 INPUT error: the value of NTI is too small; the number of necessary output points is greater than NTI-1. TERMINAL ERROR. 123 INPUT error: the value of LWG is less than the number of output points. Increase the dimension of the array WG and the value of LWG. TERMINAL ERROR. 216 This indicates that during integration the requested accuracy could not be achieved. User must increase error tolerance. TERMINAL ERROR. 219 This indicates that the routine needs more space to store the integra- tion grid point. An estimate for the required workspace (i.e. the value for LWG) is given. TERMINAL ERROR. 230 This indicates that Newton failed to converge. TERMINAL ERROR. 231 This indicates that the number of iteration has become greater than ITLIM. TERMINAL ERROR. 240 This indicates that the global error is probably larger than the error tolerance due to instabilities in the system. Most likely the system is ill-conditioned. Output value is the estimated error amplification factor. WARNING ERROR. 250 This indicate that one of the upper triangular matrices U is singular. TERMINAL ERROR. 260 This indicates that the problem is probably too ill-conditioned with respect to the BC. TERMINAL ERROR. **************** Auxiliary routines **************** Calls are made to the MUS library routines: DBCMAV, DCHINC, DCROUT, DCSAOJ, DCSHPO, DFQUS, DFUNRC, DGTUR, DINPRO, DINTCH, DJINGX, DKPCH, DLUDEC, DMATVC, DNEWPO, DPSR, DQEVAL, DQUDEC, DRKF1S, DRKFGG, DRKFGS, DRKFMS, DSBVP, DSOLDE, DSOLUP, DSORTD, DTAMVC, ERRHAN. **************** Remarks **************** MUSN is written by R.M.M. Mattheij and G.W.M. Staarink. Last update: 02-15-88. **************** Example of the use of MUSN **************** Consider the differential equation: u' = .5 u*(w-u) / v v' = -0.5 (w-u) w' = (0.9 - 1000 (w-y) - 0.5 w(w-u)) / x x' = 0.5 (w-u) y' = -100 (y-w) and the boundary conditions: u(0) = v(0) = w(0) = 1 x(0) = -10 w(1) = y(1) As an initial guess for the solution we take: u(t) = 1 ; v(t) = 1 ; x(t) = -10 ; w(t) = -4.5 t*t + 8.91 t + 1 ; y(t) = -4.5 t*t + 9 t + 0.91 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ER(5),TI(12),X(5,12),Q(5,5,12),U(15,12),D(5,12), 1 PHIREC(15,12),W(315),WGR(20) INTEGER IW(27) EXTERNAL FDIF,Y0T,G C C SETTING OF THE INPUT PARAMETERS C N = 5 NU = 15 NTI = 12 LW = 315 LIW = 27 ER(3) = 1.1D-15 LWG = 20 ER(1) = 1.D-6 ER(2) = 1.D-2 A = 0.D0 B = 1.D0 NRTI = 10 AMP = 100 ITLIM = 20 CALL MUSN(FDIF,Y0T,G,N,A,B,ER,TI,NTI,NRTI,AMP,ITLIM,X,Q,U,NU,D, 1 PHIREC,KPART,W,LW,IW,LIW,WGR,LWG,IERROR) . . . END SUBROUTINE FDIF(T,Y,F) C ---------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(5),F(5) C Y3MY1 = Y(3) - Y(1) Y3MY5 = Y(3) - Y(5) F(1) = 0.5D0 * Y(1) * Y3MY1 / Y(2) F(2) = - 0.5D0 * Y3MY1 F(3) = (0.9D0 - 1.D3 * Y3MY5 - 0.5D0 * Y(3) * Y3MY1) / Y(4) F(4) = 0.5D0 * Y3MY1 F(5) = 1.D2 * Y3MY5 RETURN C END OF FDIF END SUBROUTINE Y0T(T,X) C ------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(5) C X(1) = 1.D0 X(2) = 1.D0 X(4) = -10.D0 X(3) = - 4.5D0*T*T + 8.91D0 * T + 1.D0 X(5) = - 4.5D0*T*T + 9.D0 * T + 0.91D0 RETURN C END OF Y0T END SUBROUTINE G(N,XA,XB,FG,DGA,DGB) C -------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XA(N),XB(N),FG(N),DGA(N,N),DGB(N,N) C DO 1100 I = 1 , N DO 1100 J = 1 , N DGA(I,J) = 0.D0 DGB(I,J) = 0.D0 1100 CONTINUE DGA(1,1) = 1.D0 DGA(2,2) = 1.D0 DGA(3,3) = 1.D0 DGA(4,4) = 1.D0 DGB(5,3) = 1.D0 DGB(5,5) = -1.D0 FG(1) = XA(1) - 1.D0 FG(2) = XA(2) - 1.D0 FG(3) = XA(3) - 1.D0 FG(4) = XA(4) + 10.D0 FG(5) = XB(3) - XB(5) RETURN C END OF G END