Main integration routine for the rk_variable_step_class.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(rk_variable_step_class), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in) | :: | t0 |
initial time |
||
| real(kind=wp), | intent(in), | dimension(:) | :: | x0 |
initial state |
|
| real(kind=wp), | intent(in) | :: | h |
initial abs(time step) |
||
| real(kind=wp), | intent(in) | :: | tf |
final time |
||
| real(kind=wp), | intent(out), | dimension(:) | :: | xf |
final state |
|
| integer, | intent(out), | optional | :: | ierr |
0 = no errors, <0 = error. if not present, an error will stop program. |
subroutine integrate(me,t0,x0,h,tf,xf,ierr) implicit none class(rk_variable_step_class),intent(inout) :: me real(wp),intent(in) :: t0 !! initial time real(wp),dimension(:),intent(in) :: x0 !! initial state real(wp),intent(in) :: h !! initial abs(time step) real(wp),intent(in) :: tf !! final time real(wp),dimension(:),intent(out) :: xf !! final state integer,intent(out),optional :: ierr !! 0 = no errors, !! <0 = error. !! if not present, an error will stop program. real(wp) :: t,dt,t2,err,tol,dt_new real(wp),dimension(me%n) :: x,terr,etol,xp0 logical :: last,export,accept integer :: i,p if (present(ierr)) ierr = 0 if (.not. associated(me%f)) then if (present(ierr)) then ierr = -1 return else error stop 'Error in integrate: f is not associated.' end if end if me%num_rejected_steps = 0 export = associated(me%report) if (export) call me%report(t0,x0) !first point if (t0==tf) then xf = x0 else t = t0 x = x0 if (h==zero) then ! compute an appropriate initial step size: ! WARNING: this may not be working in all cases ..... etol = me%rtol * me%stepsize_method%norm(x0) + me%atol call me%f(t0,x0,xp0) ! get initial dx/dt select case (me%hinit_method) case(1) call me%hstart(t0,tf,x0,xp0,etol,dt) case(2) dt = me%hinit(t0,x0,sign(one,tf-t0),xp0,me%stepsize_method%hmax,me%atol,me%rtol) case default if (present(ierr)) then ierr = -2 return else error stop 'invalid hinit_method selection' end if end select !write(*,*) 'inital step size: ',dt else ! user-specified initial step size: dt = sign(h,tf-t0) ! (correct sign) end if p = me%order() !order of the method do t2 = t + dt last = ((dt>=zero .and. t2>=tf) .or. & !adjust last time step (dt<zero .and. t2<=tf)) ! if (last) dt = tf-t ! do i=0,me%stepsize_method%max_attempts ! take a step: call me%step(t,x,dt,xf,terr) ! evaluate error and compute new step size: err = me%stepsize_method%norm(terr) tol = me%stepsize_method%norm( me%rtol * xf + me%atol ) call me%stepsize_method%compute_stepsize(dt,tol,err,p,dt_new,accept) dt = dt_new if (accept) then !accept this step exit else !step is rejected, repeat step with new dt me%num_rejected_steps = me%num_rejected_steps + 1 !note: if we have reached the min step size, and the error !is still too large, we can't proceed. if (i>=me%stepsize_method%max_attempts) then if (present(ierr)) then ierr = -3 return else error stop 'error: too many attempts to reduce step size.' end if end if if (abs(dt) <= abs(me%stepsize_method%hmin)) then if (present(ierr)) then ierr = -4 return else error stop 'warning: min step size.' end if end if !...... !... if we have two rejected steps and the step size hasn't changed.. ! then we need to abort, since no progress is being made... !...... last = ((dt>=zero .and. t2>=tf) .or. & !adjust last time step (dt<zero .and. t2<=tf)) ! if (last) dt = tf-t ! t2 = t + dt end if end do if (last) exit if (export) call me%report(t2,xf) !intermediate point x = xf t = t2 end do end if if (export) call me%report(tf,xf) !last point end subroutine integrate