! Declarations Call linear algebra setup routine Call integrator setup routine irevcm=0 1000 Call d02nnf(neq, neqmax, t, tout, y, ydot, rwork, rtol, atol, itol, inform, ysave, ny2dim, wkjac, nwkjac, jacpvt, njcpvt, imon, inln, ires, irevcm, lderiv, itask, itrace, ifail) If (irevcm.gt.0) Then If (irevcm.gt.7 .and. irevcm.lt.11) Then If (irevcm.eq.8) Then supply the Jacobian matrix (i) Else If (irevcm.eq.9) Then perform monitoring tasks requested by the user (ii) Else If (irevcm.eq.10) Then indicates an unsuccessful step End If Else evaluate the residual (iii) Endif Go To 1000 End If ! post processing (optional linear algebra diagnostic call ! (sparse case only), optional integrator diagnostic call) Stop End
|Supply the Jacobian matrix
You need only provide this facility if the argument (or if using sparse matrix linear algebra) in a call to the linear algebra setup routine (see jceval in d02nuf). If the Jacobian matrix is to be evaluated numerically by the integrator, then the remainder of section (i) can be ignored.
We must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, , has the form
where is the current step size and is an argument that depends on the integration method in use. The vector is the current solution and the vector depends on information from previous time steps. This means that .
The system of nonlinear equations that is solved has the form
but is solved in the form
where is the function defined by
It is the Jacobian matrix that you must supply as follows:
where , and are located in , and respectively and the arrays y and ydot contain the current solution and time derivatives respectively. Only the nonzero elements of the Jacobian need be set, since the locations where it is to be stored are preset to zero.
Hereafter in this document this operation will be referred to as jac.
|Perform tasks requested by you
This operation is essentially a monitoring function and additionally provides the opportunity of changing the current values of y, ydot, hnext (the step size that the integrator proposes to take on the next step), hmin (the minimum step size to be taken on the next step), and hmax (the maximum step size to be taken on the next step). The scaled local error at the end of a time step may be obtained by calling the real function d02zaf as follows:
ifail = 1 errloc = d02zaf(neq,rowk(51+neqmax),rwork(51),ifail) ! CHECK IFAIL BEFORE PROCEEDING
The following gives details of the location within the array rwork of variables that may be of interest to you:
You are advised to consult the description of monitr in d02ngf for details on what optional input can be made.
If either y or ydot are changed, then imon must be set to before return to d02nnf. If either of the values hmin or hmax are changed, then imon must be set before return to d02nnf. If hnext is changed, then imon must be set to before return to d02nnf.
In addition you can force d02nnf to evaluate the residual vector
by setting and and then returning to d02nnf; on return to this monitoring operation the residual vector will be stored in , for .
Hereafter in this document this operation will be referred to as monitr.
|Evaluate the residual
This operation must evaluate the residual
in one case and the reduced residual
in another, where is located in . The form of the residual that is returned is determined by the value of ires returned by d02nnf. If , then the residual defined by equation (2) above must be returned; if , then the residual returned by equation (1) above must be returned.
Hereafter in this document this operation will be referred to as resid.
Call d02nrf(j, iplace, inform)
Program Text (d02nnfe.f90)
Program Data (d02nnfe.d)
Program Results (d02nnfe.r)