Skip to content

Commit

Permalink
Implement UnboundedTaylorPicardIntegrator that avoids use of Bounder, r…
Browse files Browse the repository at this point in the history 
…esolve #622
  • Loading branch information
@lgeretti
lgeretti committed Jun 3, 2021
1 parent 7ab5e85 commit 9b3a855
Show file tree
Hide file tree
Showing 2 changed files with 181 additions and 8 deletions.
126 changes: 119 additions & 7 deletions source/solvers/integrator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,19 @@ IntegratorBase::flow_step(const ValidatedVectorMultivariateFunction& vf, const E
}
}

FlowStepModelType
UnboundedIntegratorBase::flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& dx, StepSizeType& hmax) const
{
StepSizeType& h=hmax;
while(true) {
try {
return flow_step(vf,dx,h,dx);
} catch(const FlowTimeStepException& e) {
h=hlf(h);
}
}
}

inline ExactDouble operator*(ExactDouble, TwoExp);
inline ExactDouble operator/(ExactDouble, TwoExp);

Expand Down Expand Up @@ -225,15 +238,9 @@ TaylorPicardIntegrator::_flow_step(const ValidatedVectorMultivariateFunction& f,
ARIADNE_THROW(FlowTimeStepException,"TaylorPicardIntegrator::flow_step","Integration of "<<f<<" starting in "<<D<<" over time interval "<<T<<" of length "<<h<<" has error "<<phi.error()<<" after "<<this->_maximum_temporal_order<<" iterations, which exceeds step maximum error "<<this->step_maximum_error());
}

FlowStepModelType res=restriction(phi,dom);

ARIADNE_LOG_PRINTLN("result="<<res);
return res;

return phi;
}



Void TaylorPicardIntegrator::_write(OutputStream& os) const {
os << "TaylorPicardIntegrator"
<< "(maximum_error = " << this->maximum_error()
Expand All @@ -246,6 +253,111 @@ Void TaylorPicardIntegrator::_write(OutputStream& os) const {
<< " )";
}

UnboundedTaylorPicardIntegrator::UnboundedTaylorPicardIntegrator(MaximumError err, Order order)
: UnboundedIntegratorBase(err,GradedSweeper<FloatDP>(DoublePrecision(),order),LipschitzConstant(0.5)),
_step_maximum_error(cast_exact(err.value())), _error_refinement_minimum_improvement_percentage(cast_exact(0.02)), _sweeper(GradedSweeper<FloatDP>(DoublePrecision(), order)), _order(order) { }

FlowStepModelType
UnboundedTaylorPicardIntegrator::flow_step(const ValidatedVectorMultivariateFunction& vf, const ExactBoxType& dx, const StepSizeType& h, const UpperBoxType& bx) const
{
ARIADNE_PRECONDITION(vf.result_size()==dx.dimension());
ARIADNE_PRECONDITION(vf.argument_size()==dx.dimension());
ARIADNE_PRECONDITION(bx.dimension()==dx.dimension());
return this->_flow_step(vf,dx,IntervalDomainType(0,h),BoxDomainType(0u),bx);
}

FlowStepModelType
UnboundedTaylorPicardIntegrator::flow_step(const ValidatedVectorMultivariateFunction& f, const ExactBoxType& D, const Interval<StepSizeType>& T, const ExactBoxType& A, const UpperBoxType& B) const
{
ARIADNE_PRECONDITION(f.result_size()==D.dimension());
ARIADNE_PRECONDITION(f.argument_size()==D.dimension()+T.dimension()+A.dimension());
ARIADNE_PRECONDITION(B.dimension()==D.dimension());
return this->_flow_step(f,D,IntervalDomainType(T),A,B);
}

FlowStepModelType
UnboundedTaylorPicardIntegrator::_flow_step(const ValidatedVectorMultivariateFunction& f, const ExactBoxType& D, const ExactIntervalType& T, const ExactBoxType& A, const UpperBoxType& B) const
{
ARIADNE_LOG_SCOPE_CREATE;
ARIADNE_LOG_PRINTLN("f="<<f);
ARIADNE_LOG_PRINTLN("D="<<D<<" T="<<T<<", A="<<A<<", B="<<B);

const bool is_autonomous = (f.argument_size()==D.dimension()+A.dimension());

const SizeType nx=D.size();
const SizeType na=A.size();

Range tarng = is_autonomous ? Range(nx+1u,nx+1u+na) : Range(nx,nx+1u+na);

StepSizeType t=static_cast<StepSizeType>(T.lower_bound());
StepSizeType h=static_cast<StepSizeType>(T.upper_bound())-t;

// Time interval centred on initial time, which will make the antiderivative more efficient
ExactIntervalType wT(t-h,t+h);
ARIADNE_ASSERT(t==med(wT));

ExactBoxType dom=join(D,T,A);
ExactBoxType wdom=join(D,wT,A);
UpperBoxType const& bx=B;
ARIADNE_LOG_PRINTLN_AT(2,"dom="<<dom<<", wdom="<<wdom);

FlowStepModelType phi0=this->function_factory().create_projection(wdom,range(0,nx));
ARIADNE_LOG_PRINTLN_AT(1,"phi0="<<phi0);
FlowStepModelType phi=this->function_factory().create_constants(wdom,cast_singleton(bx));
FlowStepModelType ta=this->function_factory().create_projection(wdom,tarng);

ARIADNE_LOG_PRINTLN_AT(1,"phi="<<phi);
for (DegreeType k=0; k!=this->_order; ++k) {
FlowStepModelType fphi=compose(f,join(std::move(phi),ta));
ARIADNE_LOG_PRINTLN_AT(2,"fphi="<<fphi);
phi=antiderivative(fphi,nx)+phi0;
ARIADNE_LOG_PRINTLN_AT(2,"phi="<<phi);
}
auto errors = phi.errors();
ARIADNE_LOG_PRINTLN_AT(2,"initial errors to validate=" << errors);
FlowStepModelType fphi=compose(f,join(std::move(phi),ta));
phi=antiderivative(fphi,nx)+phi0;
auto new_errors = phi.errors();
for (SizeType i=0; i<errors.size(); ++i) {
if (not refines(new_errors[i],errors[i])) {
ARIADNE_THROW(FlowTimeStepException,"UnboundedTaylorPicardIntegrator::flow_step","Integration of "<<f<<" starting in "<<D<<" over time interval "<<T<<" of length "<<h<<" has errors "<<new_errors<<", which are not smaller than previous errors " << errors);
}
}
errors = new_errors;
ARIADNE_LOG_PRINTLN_AT(2,"validated errors=" << errors);
while (true) {
fphi=compose(f,join(std::move(phi),ta));
phi=antiderivative(fphi,nx)+phi0;
new_errors = phi.errors();
Bool has_improved = false;
for (SizeType i=0; i<errors.size(); ++i) {
auto error_improvement = cast_exact((errors[i]-new_errors[i])/errors[i]);
if (error_improvement >= this->_error_refinement_minimum_improvement_percentage) {
has_improved = true;
break;
}
}
if (not has_improved) break;
errors = new_errors;
ARIADNE_LOG_PRINTLN_VAR_AT(2,errors);
}

if (phi.error().raw()>this->maximum_error()) {
ARIADNE_THROW(FlowTimeStepException,"UnboundedTaylorPicardIntegrator::flow_step","Integration of "<<f<<" starting in "<<D<<" over time interval "<<T<<" of length "<<h<<" has error "<<phi.error()<<", which exceeds maximum error "<<this->maximum_error());
}

return phi;
}

Void UnboundedTaylorPicardIntegrator::_write(OutputStream& os) const {
os << "UnboundedTaylorPicardIntegrator"
<< "(maximum_error = " << this->maximum_error()
<< ", function_factory = " << this->function_factory()
<< ", order = " << this->order()
<< " )";
}


} // namespace Ariadne


Expand Down
63 changes: 62 additions & 1 deletion source/solvers/integrator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,21 @@ class IntegratorBase
BounderPointer _bounder_ptr;
};

class UnboundedIntegratorBase : public IntegratorBase {
protected:
UnboundedIntegratorBase(MaximumError e, GradedSweeper<FloatDP> sweeper, LipschitzConstant lipschitz) :
IntegratorBase(e, sweeper, lipschitz) { }

public:

virtual FlowStepModelType
flow_step(const ValidatedVectorMultivariateFunction& vector_field,
const ExactBoxType& state_domain,
StepSizeType& suggested_time_step) const override;

using IntegratorBase::flow_step;
};

//! \brief An integrator which uses a validated Picard iteration on Taylor models.
class TaylorPicardIntegrator
: public IntegratorBase
Expand All @@ -189,7 +204,6 @@ class TaylorPicardIntegrator
Void set_maximum_temporal_order(Nat m) { this->_maximum_temporal_order=m; }
//! \brief Set the sweep threshold of the Taylor model.
Sweeper<FloatDP> const& sweeper() const { return this->_sweeper; }
Void set_sweeper(Sweeper<FloatDP> const& sweeper) { _sweeper = sweeper; }
//! \brief Set the maximum error of a single step.
ExactDouble step_maximum_error() const { return this->_step_maximum_error; }
Void set_step_maximum_error(ApproximateDouble e) { _step_maximum_error = cast_exact(e); }
Expand Down Expand Up @@ -221,6 +235,53 @@ class TaylorPicardIntegrator
const UpperBoxType& bounding_box) const;
};

class UnboundedTaylorPicardIntegrator
: public UnboundedIntegratorBase
{
ExactDouble _step_maximum_error;
ExactDouble _error_refinement_minimum_improvement_percentage;
GradedSweeper<FloatDP> _sweeper;
DegreeType _order;
public:
//! \brief Default constructor.
UnboundedTaylorPicardIntegrator(MaximumError err, Order order);

//! \brief The order of the method.
DegreeType order() const { return this->_order; }
Void set_order(Nat m) { this->_order=m; this->_sweeper = GradedSweeper<FloatDP>(DoublePrecision(),m); }
//! \brief Set the maximum error of a single step.
ExactDouble step_maximum_error() const { return this->_step_maximum_error; }
Void set_step_maximum_error(ApproximateDouble e) { _step_maximum_error = cast_exact(e); }
ExactDouble error_refinement_minimum_improvement_percentage() const { return this->_error_refinement_minimum_improvement_percentage; }
Void set_error_refinement_minimum_improvement_percentage(ApproximateDouble e) { _error_refinement_minimum_improvement_percentage = cast_exact(e); }

virtual UnboundedTaylorPicardIntegrator* clone() const { return new UnboundedTaylorPicardIntegrator(*this); }
virtual Void _write(OutputStream& os) const;

virtual FlowStepModelType
flow_step(const ValidatedVectorMultivariateFunction& vector_field,
const ExactBoxType& state_domain,
const StepSizeType& time_step,
const UpperBoxType& bounding_box) const;

virtual FlowStepModelType
flow_step(const ValidatedVectorMultivariateFunction& differential_equation,
const ExactBoxType& state_domain,
const Interval<StepSizeType>& time_domain,
const ExactBoxType& parameter_domain,
const UpperBoxType& bounding_box) const;

using UnboundedIntegratorBase::flow_step;

private:
FlowStepModelType
_flow_step(const ValidatedVectorMultivariateFunction& vector_field_or_differential_equation,
const ExactBoxType& state_domain,
const ExactIntervalType& time_domain,
const ExactBoxType& parameter_domain,
const UpperBoxType& bounding_box) const;
};



//! \brief An integrator which computes the Taylor series of the flow function with remainder term.
Expand Down

0 comments on commit 9b3a855

Please sign in to comment.