Name : MIRKDC Authors: Wayne Enright (enright@cs.toronto.edu), Paul Muir (muir@stmarys.ca). Additional programming assistance: Mark Adams, Nicole DeYoung. Thanks to Beta-Testers: Ian Gladwell, Larry Shampine, Richard Pancer. Note: This code makes use of the Level 1 BLAS routines. These must be linked with the MIRKDC package in order for it to work. Purpose : Solves a system of first order, non-linear, boundary value, ordinary differential equations, y'(t) = f(t,y(t)), with separated boundary conditions, g_a(y(a))=0 and g_b(y(b))=0. Based on a mesh of points which subdivides [a,b], the code uses mono-implicit Runge-Kutta (MIRK) methods to discretize the ODE's and solves the resultant non-linear algebraic equations, using a hybrid damped Newton/fixed Jacobian iteration, to get a discrete solution approximation at the meshpoints. Continuous MIRK schemes are used to augment the discrete solution with a C^1 continuous interpolant, u(t), which is then used to compute an estimate of the solution defect, | u'(t) - f(t,u(t)) |. The user defined tolerance is then applied to the defect as the termination criterion. Adaptive mesh redistribution is performed to equidistribute the estimate of the defect. For further details see [Enright,Muir, "Runge-Kutta Software with Defect Control for Boundary Value ODES",SIAM J. Sci. Comput. 17 (1996),479-497]. Usage : call mirkdc(method, tol, neqns, Nsub, MxNsub, mesh_input, mesh, sol_input, Y, ldY, output_control, info, iwork, work, init_de, init_Y, fsub, gsub, dfsub, dgsub) Arguments : method - (input) - integer. Possible values for `method' and the corresponding MIRK schemes are: Value of `method' | MIRK formula ---------------------------------------------------------- 2 | Second order MIRK scheme | (Trapezoidal scheme). 4 | Fourth order MIRK scheme | (Lobatto collocation scheme). 6 | Sixth order MIRK scheme. tol - (input) - double precision. A user-defined tolerance applied to an estimate of the defect of the approximate solution; the defect is the amount by which the continuous approximate solution fails to satisfy the ODE system. MIRKDC attempts to return, a continuous solution approximation, u(t), such that |u'(t)-f(t,u(t))|/(|f(t,u(t))|+1) is less than `tol'. (This is a blended absolute/relative measure of the defect.) The same tolerance is applied to each component of the defect. neqns - (input) - integer. The number of first order differential equations. This is also the total number of boundary conditions. Nsub - (input/output) - integer. When passed to MIRKDC it is the initial user-defined number of subintervals into which the problem interval is partitioned. Upon return from MIRKDC it is the number of subintervals in the final mesh. MxNsub - (input) - integer. The user-defined maximum allowable number of subintervals. The value of this parameter will influence storage requirements - see e.g. `iwork' and `work'. mesh_input - (input) - integer. a) A value of 0 indicates that a uniform initial mesh is to be set up by MIRKDC. b) A value of 1 indicates that an initial mesh is already available in the array 'mesh'. This option should be chosen if the user knows which regions of the interval correspond to difficult solution behavior, e.g. a boundary layer; the mesh should be finer in those regions. c) A value of 2 indicates that the user is employing a ``continuation sequence'' of problems and that the mesh from a previous problem will be available in the array `mesh'. mesh - (input/output) - double precision array `(0:MxNsub)'. On input to MIRKDC this is the set of points that partition the problem interval ( empty if `mesh_input' is 0 ). On successful return from MIRKDC it is the set of points defining the final mesh. sol_input - (input) - integer. a) A value of 1 indicates that the initial solution estimate will be provided by the user through the `init_Y' routine. This is the usual value for `sol_input' and the usual mode of execution for MIRKDC. It is recommended that `Y' not be set to a constant in case there is a singularity in f(t,y(t)), for that constant value. b) A value of 2 indicates that the initial solution estimate at each mesh point is available in the array `Y'. This option should only be chosen when solving a continuation sequence of problems, where `Y' takes its value from the solution of a previous problem. Y - (input/output) - double precision two dimensional array of size `ldY x (MxNsub+1)'. On input to MIRKDC, if `sol_input'= 2, this is the discrete approximation to the solution, evaluated at each point of the mesh contained in `mesh'. (Empty if `sol_input'=1.) The i-th column of `Y' contains the solution approximation vector for the i-th meshpoint. On successful return from MIRKDC, this is the discrete approximation to the solution evaluated at the points of the final mesh. ldY - (input) - integer. The leading dimension of the array `Y'. (`ldY' >= `neqns'). output_control - (input) - integer. Controls the amount of information output. Newton iteration counts, mesh selection information, and defect estimates are examples of such information. Possible values of `output_control' are : 0 - No output, 1 - Intermediate output, 2 - Full output. info - (output) - integer. A communication flag: `info' = 0 implies successful termination; the solution has been obtained such that the estimates of the defect satisfy the user defined tolerance - see `tol'. `info' = 1 implies unsuccessful termination; the size of the new mesh would be too large. Unless `MxNsub' is actually fairly small, this is usually an indication that the code encountered difficulty with the convergence of the Newton iterations. The user-supplied routines `fsub', `dfsub', `gsub', and `dgsub' should be checked. iwork - (output) - integer (work) array of size at least `neqns * (MxNsub + 1)'. work - (output) - double precision (work) array of size at least `MxNsub * neqns * (2*neqns+35) + 6 * neqns * (4*neqns+1)'. init_de - User-supplied subroutine to initialize the problem dependent variables, 'leftbc', 'a', and 'b'. call init_de(leftbc, a, b), where leftbc - (output) - integer. Number of boundary conditions applied at the left endpoint of the problem interval. a - (output) - double precision. The left endpoint of the problem interval. b - (output) - double precision. The right endpoint of the problem interval. init_Y - User-supplied subroutine to initialize the discrete solution approximation, `Y'. (Needed if `sol_input'=1). call init_Y(Nsub, neqns, mesh, Y, ldY), where Nsub - (input) - integer. The number of subintervals in the initial mesh. neqns - (input) - integer. The number of differential equations in the first order system. mesh - (input) - double precision array (0:Nsub). Points making up the partition of the problem interval. Y - (output) - double precision two dimensional array of size `ldY x (Nsub + 1)'. The initial approximation to the discrete solution at the points in the initial mesh. The i-th column of `Y' is the solution approximation vector at mesh point, `mesh(i-1)'. ldY - (input) - integer. The leading dimension of Y. (`ldY' >= `neqns'). fsub - User-supplied subroutine which defines f(t,y) for the first order system of differential equations, y' = f(t,y). call fsub(neqns, t, y, f), where neqns - (input) - integer. The number of differential equations in the first order system. t - (input) - double precision. A point in the problem interval. y - (input) - double precision array `(1:neqns)'. The current solution approximation at t. f - (output) - double precision array `(1:neqns)'. The value of f(t,y). gsub - User-supplied subroutine which defines the boundary condition equations. call gsub(neqns, ya, yb, bc), where neqns - (input) - integer. The number of differential equations in the first order system. Also the total number of boundary conditions. ya - (input) - double precision array `(1:neqns)'. Current solution approximation at the left endpoint. yb - (input) - double precision array `(1:neqns)'. Current solution approximation at the right endpoint. bc - (output) - double precision array `(1:neqns)'. Boundary condition equations evaluated at `ya' and `yb'. The first `leftbc' components of `bc' are `g_a(ya)'; the remaining `neqns-leftbc' components of `bc' are `g_b(yb)'. dfsub - User-supplied subroutine which defines the Jacobian, df/dy, of the system of differential equations. Since the array `JJ' has already been set to zero, it is only necessary to set the non-zero elements. call dfsub(neqns, t, y, JJ), where neqns - (input) - integer. The number of differential equations in the first order system. t - (input) - double precision. A point in the problem interval. y - (input) - double precision array `(1:neqns)'. The current solution approximation at t. JJ - (output) - double precision two dimensional array of size `neqns x neqns'. Contains the Jacobian, df/dy. dgsub - User-supplied subroutine which defines the Jacobian of the boundary conditions, that is d(bc)/dy. The array `dbc' has been set to zero prior to the call to 'dgsub', so it is only necessary for the user to set the non-zero elements. call dgsub(neqns, ya, yb, dbc), where neqns - (input) - integer. The number of differential equations in the first order system. Also the total number of boundary conditions. ya - (input) - double precision array `(1:neqns)'. Current solution approximation at the left endpoint. yb - (input) - double precision array `(1:neqns)'. Current solution approximation at the right endpoint. dbc - (output) - double precision two dimensional array of size `neqns x neqns'. Contains the Jacobian of the boundary conditions. Remarks: 1. On output, both the `iwork' and `work' arrays contain information that is required in order for the user to access the continuous solution approximation using the `sol_eval' routine. 2. The code makes use of two external packages of routines. The COLROW package [Diaz,Fairweather,Keast, "COLROW and ARECO:Fortran packages for solving certain almost block diagonal linear systems by modified row and column elimination",ACM Trans. Math. Soft. 9 (1983),376-380] is included with the MIRKDC source code. The double precision version of the level 1 BLAS routines, [Lawson,Hanson,Kincaid,Krogh,"Basic linear algebra sub-programs for Fortran usage",ACM Trans. Math. Soft. 5 (1979), 308-323], is available from the netlib software repository and must be linked with the MIRKDC package. 3. The usual technique for controlling the accuracy or quality of an appropriate solution is to estimate and control the global error. MIRKDC, in contrast, estimates the defect, u'(t) - f(t,u(t)) (where u(t) is the continuous approximate solution) and applies the user defined tolerance to control the defect. See `tol' parameter documentation and [Enright,Muir, "Runge-Kutta Software with Defect Control for Boundary Value ODES",SIAM J. Sci. Comput. 17 (1996),479-497]. 4. When using continuation, it should be noted that the continuation sequence must be chosen so that the Newton iterations in MIRKDC do not fail since the absence of an appropriate `init_Y' routine will be problematic. (See e.g. [Ascher, Christiansen, Russell, "Collocation software for boundary value ODEs", Math. Comp., 33 (1979), 659-679], for an example of how a BVODE code is used in a parameter continuation sequence.) For further details see also the internal documentation for MIRKDC. ------------------------------------------------------------------------------- Name : Sol_eval Purpose : This routine provides the solution approximation and its derivative at a given point within the problem interval. Usage : call sol_eval(t,z,z_prime,iwork,work) Arguments : t - (input) - double precision. The evaluation point within the problem interval. z - (output) - double precision array `(1:neqns)'. The value of the interpolant at `t'. z_prime - (output) - double precision array `(1:neqns)'. The value of the derivative of the interpolant at `t'. iwork - (input) - integer array of size at least `neqns * (MxNsub + 1)'. This array contains integer information employed or set in MIRKDC and needed by `sol_eval'. work - (input) double precision array of size at least `MxNsub * neqns * (2*neqns+35) + 6 * neqns * (4*neqns+1)'. This array contains double precision information employed or set in MIRKDC and needed by `sol_eval'. ------------------------------------------------------------------------------ Sample Driver program: c============================================================================== c ***Model of main program to use MIRKDC, for Swirling Flow III test problem*** c "Swirling Flow III", [Ascher, Mattheij, and Russell, c "Numerical Solution of Boundary Value Problems for Ordinary c Differential Equations", Classics in Applied Mathematics Series, c Society for Industrial and Applied Mathematics, Philadelphia, 1995]. c============================================================================== c c***declaration of constants*** c integer neqns, MxNsub, ldY parameter(neqns=6,MxNsub=3000,ldY=10) c c neqns - The number of differential equations. c MxNsub - The maximum number of subintervals. c ldY - The leading dimension of the matrix Y. c c***declaration of remaining parameters to mirkdc*** c integer method double precision tol integer Nsub integer mesh_input double precision mesh(0:MxNsub) integer sol_input double precision Y(ldY,(MxNsub+1)) integer output_control integer info double precision work(MxNsub*neqns*(2*neqns+35) + +6*neqns*(4*neqns+1)) c Make sure 'work' is at least c MxNsub*neqns*(2*neqns+35) + 6*neqns*(4*neqns+1) c integer iwork(neqns*(MxNsub+1)) c Make sure 'iwork' is at least neqns*(MxNsub+1) c double precision epsilon common /ode/ epsilon c Problem dependent parameter c c Value of 'method' | MIRK formula c ------------------------------------- c 2 | Second order MIRK scheme c 4 | Fourth order MIRK scheme c 6 | Sixth order MIRK scheme c c tol - A user-defined tolerance applied to the defect of the c computed solution; the defect is the amount by which c the continuous approximate solution fails to satisfy c the ODE system. If 'mirkdc' returns successfully, c |def(t)|/(|f(t,y(t))|+1) will be less than 'tol', where c y'(t)=f(t,y(t)) is the ODE and def(t)= y'(t)-f(t,y(t)). c The defect has one component per solution component. c The same tolerance is applied to each component of the c defect. c c Nsub - On input to 'mirkdc', the number of subintervals into which c the problem interval is partitioned. c On successful return from 'mirkdc', the number of sub - c intervals in the final mesh. c c mesh_input - a) A value of 0 indicates that a uniform initial c mesh is to be set up by 'mirkdc'. c c b) A value of 1 indicates that an initial mesh is c already available in the array 'mesh'. This option c should be chosen if the user knows which regions of c the interval are problematic so a finer mesh may be c used on these trouble areas. c c c) A value of 2 indicates that the user is doing a c continuation problem in which case the mesh from c a previous run is stored in the array 'mesh'. c c mesh - On input to 'mirkdc', the set of points that partition the c problem interval, initially empty if 'mesh_input' = 0 c On successful return from 'mirkdc', the set of points c defining the final mesh. c c sol_input - a) A value of 1 indicates that an initial solution c estimate is to be provided by the user through c the 'init_Y' routine. It is recommended that Y not c be set to a constant as there could be a singularity c at that point in the differential equation. c c b) A value of 2 indicates that an initial solution c estimate at each mesh point is available in the c matrix 'Y'. This option should only be chosen when c doing continuation problems where Y takes its c value from the solution of a previous run. c c Y - On input to 'mirkdc' if 'sol_input' = 2, Y is a matrix c containing the discrete approximation to the solution, c evaluated at each point of the mesh contained in 'mesh'. c On successful return from 'mirkdc', Y is a matrix c containing the discrete approximation to the solution c evaluated at the points in 'mesh'. c c c output_control - Controls the amount of information output. c Profiling of Newton iteration counts, mesh selection, c relative defect estimate are all examples of such c information. c Possible values of output_control are: c 0 - No output. c 1 - Intermediate output. c 2 - Full output. c c info - Communication flag: c c info = 0 successful termination - solution obtained c to within user specified tolerance. c c info = -1 unsuccessful termination - the size c of the new mesh would be too large. c external init_de,init_Y,fsub,gsub,dfsub,dgsub c c c See descriptions in the sample routines for further details. c c init_de - User - supplied subroutine to initialize the problem c dependent variables, 'leftbc', 'a', and 'b'. c As well this is a good place to set up any common blocks c needed to pass ODE parameter related info to any other c user defined routines. c call init_de(leftbc, a, b) where c leftbc - Number of boundary conditions supplied at the c left endpoint of the problem interval.(output) c a - The left endpoint of the interval.(output) c b - The right endpoint of the interval.(output) c c init_Y - User - supplied subroutine to initialize the discrete c solution approximation, Y. Y is a matrix of size c (neqns) x (Nsub + 1). The i'th column of Y is the c solution approximation vector for the mesh point, c mesh(i-1). c call init_Y(Nsub, neqns, mesh, Y, ldY) where c Nsub - The number of subintervals in the mesh.(input) c neqns - The number of differential equations in the c first order system.(input) c mesh - Points making up the partition of the problem c interval.(input) c Y - The initial approximation of the discrete c solution at the points in the initial mesh. c (output). c ldY - The leading dimension of Y.(input) c c fsub - User - supplied subroutine which defines f(t,y) for the c first order system of differential equations, y' = f(t,y). c call fsub(neqns, t, y, f) where c neqns - The number of differential equations in the c first order system.(input) c t - A point in the problem interval.(input) c y - The solution approximation at time, t.(input) c f - The value of f(t,y).(output) c c gsub - User - supplied subroutine which defines the boundary c condition equations. c call gsub(neqns, ya, yb, bc) where c neqns - The number of differential equations in the c first order system.(input) c ya - Denotes the boundary conditions that have c been specified at the left endpoint.(input) c yb - Denotes the boundary condidtions that have c been specified at the right endpoint.(input) c bc - The value of the boundary condition equation. c (output) c c dfsub - User - supplied subroutine which defines the Jacobian of c the system of differential equations with respect to y, c i.e. df/dy. The Jacobian has already been set to zero so c it is only necessary that the user input the non-zero c elements. c call dfsub(neqns, t, y, JJ) where c neqns - The number of differential equations in the c first order system.(input) c t - A point in the problem interval.(input) c y - The solution approximation at time, t.(input) c JJ - A two-dimensional array containing the c Jacobian.(output) c c dgsub - User - supplied subroutine which defines the Jacobian of c the boundary conditions, that is d bc/dy. The Jacobian c has been set to zero prior to the call to 'dgsub', so it c is only necessary that the user input the non-zero c elements. c call dgsub(neqns, ya, yb, dbc) where c neqns - The number of differential equations in the c first order system.(input) c ya - Denotes the boundary conditions that have been c specified at the left endpoint.(input) c yb - Denotes the boundary conditions that have been c specified at the right endpoint.(input) c dbc - A two-dimensional array containing the Jacobian c of the boundary conditions.(output) c c ***declaration of local variables*** c integer i,j,fail double precision x,h double precision z(neqns),zp(neqns) c c i,j - loop indexes c c fail - flag to control execution of 'sol_eval' c c x - local variable to store a meshpoint c c h - local variable to store the size of a subinterval c c z - value of the interpolant evaluated at a mesh point c c zp - value of the interpolant's derivative evaluated at a c mesh point c------------------------------------------------------------------------- c***Setup parameters*** c write(*,*) write(*,*)'Initialization :' write(*,*) c c Set the order of the MIRK scheme. write(*,*) write(*,*)'Input value of method 2, 4 or 6' read(5,*) method c c Set the defect tolerance. write(*,*)'Input value of tolerance' read(5,*)tol c c Set the number of subintervals in the initial mesh. write(6,*)'Input the initial number of subintervals.' Read(5,*) Nsub c c Set the value of problem dependent parameter, epsilon write(6,*)'Input the value of epsilon' read(5,*)epsilon c c Set mesh_input to indicate that mirkdc should begin with a uniform mesh. mesh_input = 0 c c Set indicator for type of initial solution estimate. sol_input = 1 c c Set value for output_control output_control = 1 c write(6,*) write(6,*)'Swirling Flow III' write(6,*) c c ******Call mirkdc****** c call mirkdc(method,tol,neqns,Nsub,MxNsub,mesh_input, + mesh,sol_input,Y,ldY,output_control,info,iwork,work, + init_de,init_Y,fsub,gsub,dfsub,dgsub) c c ******Return from mirkdc****** c c Upon successful return from mirkdc, print out solution. c Solution interpolant is available through the sol_eval routine c This routine requires access to the work array. c c Check return from MIRKDC c if (info .NE. 0) then write(6,*)'Unsuccessful termination' stop endif c c Compute solution approximations at 11 points. c h = (mesh(Nsub) - mesh(0))/10.0d0 do 40 i = 1, 11 x = (mesh(0) + (i-1)*h) call sol_eval(x,z,zp,fail,iwork,work) if (fail.EQ.0) then write(6,*)'At meshpoint x=', x,' the solution is' do 35 j=1,neqns write(6,100)z(j) 100 format(1x,d30.15) 35 continue else write(6,*)'Attempted to evaluate the interpolant at a point', + 'that was outside of the problem interval. No information', + 'is available for that point.' end if 40 continue stop end c A Test Problem for Mirkdc c============================================================================== c CONTENTS c c This file contains the routines associated with the nonlinear c problem "Swirling Flow III", [Ascher, Mattheij, and Russell, c "Numerical Solution of Boundary Value Problems for Ordinary c Differential Equations", Classics in Applied Mathematics Series, c Society for Industrial and Applied Mathematics, Philadelphia, 1995]. c This is the special case of counter rotating disks - thus Sigma_0=-1 c and Sigma_1=1. It is stated as: c c Original Form First Order System Form c ------------- ----------------------- c y1' = y2, c g'' = gf' - fg' c --------- , y2' = y1*y4 - y3*y2 c epsilon ------------- , c epsilon c f'''' = -ff''' - gg' y3' = y4 , c ------------ , c epsilon y4' = y5 , y5' = y6 , c c y6' = -y3*y6 - y1*y2 c -------------- , c epsilon c with c c g(0) = -1, g(1) = 1, y1(0) = -1, y1(1) = 1, c c f(0) = f'(0) = f(1) = f'(1) = 0, y3(0) = y4(0) = 0, c y3(1) = y4(1) = 0, c and c epsilon = {a small parameter}. c c No closed form solution is available. c============================================================================== c subroutine init_de(leftbc,a,b) c The purpose of this routine is to initialize the problem dependent c variables 'leftbc', the number of boundary conditions at the left c endpoint of the problem interval, 'neqns', the number of differential c equations, 'a', the left endpoint of the problem interval, and 'b', c the right endpoint of the problem interval. c------------------------------------------------------------------------------ c c***declaration of parameters*** c exports: integer leftbc double precision a,b c 'leftbc' the number of boundary conditions imposed at the c left end of the problem interval c 'a' the left endpoint of the problem interval c 'b' the right endpoint of the problem intervalc c------------------------------------------------------------------------------ c leftbc = 3 c a = 0.0d0 b = 1.0d0 c return end c c============================================================================= c subroutine init_Y(Nsub,neqns,mesh,Y,ldY) c c The purpose of this routine is to initialize the discrete solution c approximation, 'Y'. The input, 'Nsub', is the number of subintervals c in the initial partition of the problem interval. The points of the c partition are stored in 'mesh'. 'Y' has 'neqns' components per meshpoint c----------------------------------------------------------------------------- c c***declaration of parameters*** c imports: integer Nsub, neqns, ldY double precision mesh(0:Nsub) c c 'Nsub' number of subintervals in 'mesh' c 'neqns' number of differential equations c 'mesh' points making up partition of problem interval c 'ldY' the leading dimension of Y. c c exports: double precision Y(ldY,(Nsub+1)) c c 'Y' the initial approximation to the discrete solution c at the points of the initial mesh. c c***declaration of local variables*** c integer i c c 'i' loop index from 0 to Nsub c------------------------------------------------------------------------------ c The initial approximation is based on straight lines joining the c boundary conditions for each differential equation in the first c order system. Boundary conditions are given for the first, c third and fourth equations, defining the lines y = 2t-1 (or mesh(i)) c y = 0 and y = 0. Since there are no boundary conditions for the second, c fifth, or sixth equations, the line drawn for y2 is y = 1 (the c derivative of y = t since y1' = y2 ) and y5 and y6 are initialized c by the line y = 0. c---------------------------------------------------------------------------- c do 10 i = 0, Nsub Y(1, i+1) = 2.0d0*mesh(i)-1.0d0 Y(2, i+1) = 2.d0 Y(3, i+1) = 0.d0 Y(4, i+1) = 0.d0 Y(5, i+1) = 0.d0 Y(6, i+1) = 0.d0 10 continue c return end c c============================================================================= c subroutine fsub(neqns,t,y,f) c c This routine defines f(t,y) for the first order system of c differential equations.`neqns' is the size of the system. c----------------------------------------------------------------------------- c c***declaration of parameters*** c imports: integer neqns double precision t, y(neqns) c c `neqns' Number of differential equations. c `t' A point in the problem interval where f(t,y(t)) is to be c evaluated. c `y' Current solution approximation at `t', used in the c evaluation of f(t,y(t)). c exports: double precision f(neqns) c c `f' Value of f(t,y(t)) for given `t' and `y'. c c***declaration for parameter 'epsilon'*** double precision epsilon common /ode/ epsilon c------------------------------------------------------------------------------ c f(1) = y(2) f(2) = (y(1)*y(4) - y(3)*y(2))/epsilon f(3) = y(4) f(4) = y(5) f(5) = y(6) f(6) = (-y(3)*y(6) - y(1)*y(2))/epsilon c return end c c============================================================================== c subroutine dfsub(neqns,t,y,JJ) c c This routine defines the Jacobian, d f(t,y(t))/dy, of the system of c differential equations. Since `JJ' has already been set to zero, it is c only necessary to set the non-zero elements. c------------------------------------------------------------------------------ c c***declaration of parameters*** c imports: integer neqns double precision t, y(neqns) c c `neqns' Number of differential equations. c `t' A point in the problem interval where d f(t,y(t))/dy is to c be evaluated. c `y' Current solution approximation at `t', used in the c evaluation of d f(t,y(t))/dy. c exports: double precision JJ(neqns,neqns) c c `JJ' Value of d f(t,y(t))/dy for given `t' and `y'. c c declaration for parameter 'epsilon' double precision epsilon common /ode/ epsilon c------------------------------------------------------------------------------ c JJ(1,2) = 1.0d0 c JJ(2,1) = y(4)/epsilon JJ(2,2) = -y(3)/epsilon JJ(2,3) = -y(2)/epsilon JJ(2,4) = y(1)/epsilon c JJ(3,4) = 1.0d0 c JJ(4,5) = 1.0d0 c JJ(5,6) = 1.0d0 c JJ(6,1) = -y(2)/epsilon JJ(6,2) = -y(1)/epsilon JJ(6,3) = -y(6)/epsilon JJ(6,6) = -y(3)/epsilon c return end c c============================================================================== c subroutine gsub(neqns,ya,yb,bc) c c This routine computes the 'neqns' boundary conditions, assumed c to be separated, applied at either the left endpoint, a, or the c right endpoint, b. Each boundary condition is assigned to the c corresponding component of 'bc' as follows. The code will attempt c to find 'ya' and 'yb' values (i.e. ya = y(a) and yb = y(b)) in order c to make each component of 'bc' equal to zero. This each boundary c condition must be written as some quantity set equal to 0. Thus, c the boundary condition, ya(2) = 3, would be assigned to say the c first boundary condition as, bc(1) = ya(2) - 3. c------------------------------------------------------------------------------ c c***declaration of parameters*** c imports: integer neqns double precision ya(neqns),yb(neqns) c c `neqns' Number of differential equations. c `ya' Current solution approximation at `a', used in the c evaluation of g_a(y(a)). c `yb' Current solution approximation at `b', used in the c evaluation of g_b(y(b)). c exports: double precision bc(neqns) c c `bc' The first `leftbc' locations contain g_a(y(a)). c The remaining `neqns-leftbc' locations contain g_b(y(b)). c------------------------------------------------------------------------------ c bc(1) = ya(1) + 1.0d0 bc(2) = ya(3) bc(3) = ya(4) c bc(4) = yb(1) - 1.0d0 bc(5) = yb(3) bc(6) = yb(4) c return end c c============================================================================== c subroutine dgsub(neqns,ya,yb,dbc) c c This routine defines the Jacobian, d bc/dy, of the boundary conditions. c `dbc' has already been set to zero so it is only necessary to set c the non-zero elements. c------------------------------------------------------------------------------ c c***declaration of parameters*** c imports: integer neqns double precision ya(neqns),yb(neqns) c c `neqns' Number of differential equations. c `ya' Current solution approximation at `a', used in the c evaluation of d g_a(y(a))/d ya. c `yb' Current solution approximation at `b', used in the c evaluation of d g_b(y(b))/ d yb. c exports: double precision dbc(neqns,neqns) c c `dbc' Value of d bc/dy for given `t',`ya', and `yb'. c--------------------------------------------------------------------------- c dbc(1,1) = 1.0d0 dbc(2,3) = 1.0d0 dbc(3,4) = 1.0d0 dbc(4,1) = 1.0d0 dbc(5,3) = 1.0d0 dbc(6,4) = 1.0d0 c return end