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