November 8, 2000 Fortran 95 package written by W. C. Rheinboldt Dept. of Mathematics, University of Pittsburgh Pittsburgh, PA 15260 wcrhein+@pitt.edu for the numerical solution of the following six types of quasilinear DAEs DAEQ1 : quasilinear DAE of index 1 a(t,u)u' = h(t,u) f(t,u) = 0 u(t0) = u0, such that f(t0,u0) = 0 dim(rge a)= dim(rge h) = nd, dim(rge f) = na, nu = nd na, nw = 0, dim(u) = nu, DAEN1 : nonlinear implicit DAE of index one f(t,u,u',w) = 0, u(t0) = u0, v(t0) = up0, w(t0) = w0, such that f(t0,u0,up0,w0) = 0 dim(u) = nu, dim(w) = nw, nd = nu, na = nu + nw, dim(rge f)= na The solver works internally with the DAE in the quasilinear form u' = v, f(t,u,v,w) = 0 DAEQ2 : quasilinear DAE of index 2 a(t,u)u' + b(t,u)w = h(t,u) f(t,u) = 0 u(t0) = u0, such that f(t0,u0) = 0 dim(rge a)= dim(rge h) = nd, dim(rge f)= na, dim(w) = nw, dim(u) = nu, nu + nw = nd + na NOHOL : quasilinear DAE of index 2 arising as models of mechanical multibody systems with nonholonomical contraints a(u) u" + d_{u'}f(t,u,u')^T w = h(t,u,u') f(t,u,u') = 0 u(t0) = u0, u'(t0) = up0 such that f(t0,u0,up0) = 0 dim(rge a)= dim(rge h) = nd, dim(rge f)= na, dim(w) = na, dim(u) = nd DAEQ3 : quasilinear DAE of index 3 a(t,u,u')u" + b(t,u,u')w = h(t,u,u') f(t,u) = 0 u(t0)=u0, u'(t0)=up0, such that f(t0,u0) = 0 dim(rge f)= na, dim(w) = nw, dim(rge a)= dim(rge h) = nd, dim(u) = nd+na-nw EULAG : Euler-Lagrange equations of index 3 a(u) u" + df(u)^T w = h(u,u') f(u) = 0 u(t0) = u0, such that f(u0) = 0 dim(rge a) = dim(rge h) = nd, dim(rge f) = na dim(w) = na, dim(u) = nd The solution procedure follows the algorithms developed originally in Rheinboldt, W. C. MANPACK: A set of algorithms for computations on implicitly defined manifolds. Comp. Math. w. Appl. 27, (1996), 15--28. Rheinboldt, W. C. Solving algebraically explicit DAEs with the MANPACK manifold algorithms. Comp. Math. w. Appl. 27, (1997), 31--43. For further references see also Rabier, P. J. and Rheinboldt, W. C. Nonholonomic Motion of Rigid Mechanical Systems from a DAE Viewpoint. SIAM Publications, Philadelphia, PA 2000 Rabier, P. J. and Rheinboldt, W. C. Theoretical and Numerical Analysis of Differential-Algebraic Equations Handbook of Numerical Analysis, Elsevier, The Netherlands, in print The distribution consists of the following tar-gzip files A. dae_solve.tgz containing the solver and associated libraries B. sampl_prob.tgz containing sample programs for the above six types of DAEs List of contents: A. dae_solve.tgz A.1 dae_lib.f95 This file contain the three library modules: module def_kinds to define kind-parameters for integer and real arithmetic module file_sup programs for handling file support module lin_lib a small linear algebra library A.2 dae_solver.f95 This files contains module dae_dat to set global data for the dae solver module dae_wrk to define work arrays and to provide associated routines for setting and deleting these arrays module dae_man to define local parametrizations on the manifolds and to provide the associated routines for establishing bases and for the evauation of points on the manifolds module dae_loc to provide subroutines for evaluating the local odes for the various dae types subroutine daeslv global routine for solving the dae B. sampl_prob.tgz This file contains sample programs with expected output for the following six problems catal.f95 uses DAEQ1 for a catalytic converter problem eight.f95 uses DAEN1 for a "figure eight" problem ap_q2.f95 uses DAEQ2 for quasilinear DAE given by Ascher and Petzold in 1991 tskate.f95 uses NOHOL for a DAE modelling a scating body on a tilted plane ap_q3.f95 uses DAEQ3 for a modified version of a quasilinear DAE given by Ascher and Petzold in 1991 sevbd.f95 uses EULAG for the seven body Andrews squeezing mechanism Template for the use of the package: 1. Write a module for the global data of the problem module prob_dat ! use def_kinds, ONLY : RP => RPX implicit none ! real(kind=RP) :: ...................... ! end module prob_dat 2. Write a main program for the problem program .... ! use def_kinds, ONLY : RP => RPX use dae_dat use prob_dat implicit none intrinsic ... ! interface subroutine daeslv(info) use def_kinds, ONLY : RP => RPX implicit none integer, intent(out) :: info end subroutine daeslv end interface ! ! provide the type of dae, the problem name and the ! dimensions of the problem ! integer :: info ! return indicator character(len=5) :: dnm = '.....' ! dae type character(len=6) :: pname = '......' ! problem name integer :: outunit = 10 ! output unit integer :: nvd = , & ! no. of diff. equ. nva = , & ! no. of alg. equ. nvu = , & ! dim. of u nvw = ! dim of w ! call setup(dnm,pname,outunit,nvd,nva,nvu,nvw,info) if(info /= 0) then if(info == -2) write(6,*)'unknown type of DAE' if(info == -1) write(6,*)'dimension error' write(6,*)'Stop due to input error' stop endif ! ! set initial point ! x(0:kx) = 0.0_RP .................. ! set the components of x ! ! set global data (the solver has defaults for all except tout) ! reltol = ......... ! relative tolerance abstol = ......... ! absolute tolerance h = .......... ! first step hint = .......... ! interpolation step tout = .......... ! terminal time hmax = .......... ! maximum step ! ! print any desired information about the run, such as ! write(outunit,20)reltol,abstol 20 format(' tolerances: reltol=',e14.6,' abstol= ',e14.6) write(outunit,30)h,tout 30 format(' h= ',e14.6,' tout= ',e14.6) ! call daeslv(info) ! call the solver ! call close_run(info) ! terminate run ! write(6,*)'Stop with info= ',info ! write stop message ! end program .... 3. Write a routine for the evaluation of the coefficient functions of the DAE subroutine daefct(fname,xc,vc,vec,mat,iret) ! use def_kinds, ONLY : RP => RPX use dae_dat use prob_dat implicit none intrinsic .... ! character(len=4), intent(in) :: fname real(kind=RP), dimension(0:), intent(in) :: xc real(kind=RP), dimension(0:), intent(in) :: vc real(kind=RP), dimension(0:), intent(out) :: vec real(kind=RP), dimension(0:,0:), intent(out) :: mat integer, intent(out) :: iret optional vc,vec,mat ! real(kind=RP) :: ....... ! any local variables !------------------------------------------------------ ! Routine for evaluating the parts of the dae under ! control of the identifier fname. ! Return indicator: ! iret = 0 no error ! iret = 1 not available part as requested by fname ! iret = 2 wrong dimension for output !------------------------------------------------------ ! iret = 0 select case(fname) ! case ('amt ') if(.not.present(mat)) then iret = 2 return endif mat(0,0) = .... ............. ! case ('bmt ') if(.not.present(mat)) then iret = 2 return endif mat(0,0) = .... .............. ! case ('hf ') if(.not.present(vec)) then iret = 2 return endif vec(0) = ... .............. ! case ('ff ') if(.not.present(vec)) then iret = 2 return endif vec(0) = ... .............. ! case ('dff ') if(.not.present(mat)) then iret = 2 return endif mat(0,0) = .... ............... ! case ('d2ff') if(.not.present(vc) .or. .not.present(vec)) then iret = 2 return endif vec(0) = ..... ............. ! case default iret = 1 ! end select ! end subroutine daefct The cases to be included in this subroutine differ for the six types of DAEs: DAEQ1 requires subroutine daefct with fname = 'amt ' for evaluating a(t,u) fname = 'hf ' for evaluating h(t,u) fname = 'ff ' for evaluating f(t,u) fname = 'dff ' for evaluating df(t,u) DAEN1 requires subroutine daefct with fname = 'ff ' for evaluating f(t,u) fname = 'dff ' for evaluating df(t,u) DAEQ2 requires subroutine daefct with fname = 'amt ' for evaluating a(t,u) fname = 'bmt ' for evaluating b(t,u) fname = 'hf ' for evaluating h(t,u) fname = 'ff ' for evaluating f(t,u) fname = 'dff ' for evaluating df(t,u) NOHOL requires subroutine daefct with fname = 'amt ' for evaluating a(u) fname = 'hf ' for evaluating h(t,u,u') fname = 'ff ' for evaluating f(t,u,u') fname = 'dff ' for evaluating df(t,u,u') DAEQ3 requires subroutine daefct with fname = 'amt ' for evaluating a(t,u,u') fname = 'bmt ' for evaluating b(t,u,u') fname = 'hf ' for evaluating h(t,u,u') fname = 'ff ' for evaluating f(t,u) fname = 'dff ' for evaluating df(t,u) fname = 'd2ff' for evaluating d2f(t,u)(v,v) EULAG requires subroutine daefct with fname = 'amt ' for evaluating a(u) fname = 'hf ' for evaluating h(u,u') fname = 'ff ' for evaluating f(u) fname = 'dff ' for evaluating df(u) fname = 'd2ff' for evaluating d2f(u)(v,v) 4. Write a subroutine for the intermediate output subroutine solout(nstp,xc,hinch,info) ! use def_kinds, ONLY : RP => RPX use dae_dat use prob_dat implicit none intrinsic .... ! integer, intent(in) :: nstp real(kind=RP), dimension(0:), intent(in) :: xc real(kind=RP), intent(out) :: hinch integer, intent(out) :: info ! !------------------------------------------------------------ ! Routine for intermediate printout of the solution vector xc ! ! nstp current count of accepted steps ! hinch a new value of the interpolation step or ! zero if no change is desired ! info nonzero if integration is to stop or ! zero otherwise !----------------------------------------------------------- ! ! sample output ! write(outlun,10) nstp,xc(0) 10 format(1X/' nstp = ',I6,' t= ',e18.11) ! write(outlun,20)xc(1:kc) 20 format(' xc = '/(4e16.8)) ! end subroutine solout