diff --git a/src/integrale.cpp b/src/integrale.cpp deleted file mode 100644 index c354e7ec43cb76031c972baafe309f1320135354..0000000000000000000000000000000000000000 --- a/src/integrale.cpp +++ /dev/null @@ -1,289 +0,0 @@ -#include "integrale.h" - -integrale::integrale(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function* f) : _inf_bound(inf), _sup_bound(sup), _f(f) -{ - _zeros_inf = nullptr; - _zeros_sup = nullptr; - _result = nullptr; -} - -integrale::~integrale() -{ - if (_zeros_inf) - delete _zeros_inf; - if (_zeros_sup) - delete _zeros_sup; - if (_result) - delete _result; - -// if (_f) -// delete _f; -} - - - std::list * integrale::construct_zeros_intern(const ibex::Interval & inf_bound) -{ -#ifdef __DEBUG__ - __SHOW_FCT_NAME__ -#endif - //ATTENTION pas de cout << elem << endl;preuve d'unicité - std::list * res = new std::list(); - - ibex::CtcNewton newton(*_f); - - ibex::CtcFwdBwd fb (*_f); - ibex::RoundRobin bisector(0); - - ibex::CellStack buff; - - ibex::CtcCompo comp( fb, newton); - - ibex::CtcFixPoint fix(comp, 1e-7); - - ibex::Solver s(fix, bisector, buff); - - - ibex::IntervalVector inf(1); - inf[0] = inf_bound; - std::vector sols = s.solve(inf); - - unsigned int i; - for (i = 0; i < sols.size(); i++){ - //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 - ibex::IntervalVector temp(sols[i].diam()); - if (i < sols.size() - 1){ - for (int j = 0; j < sols[i].size(); j++){ - temp[j] = (sols[i][j] + sols[i+1][j]) / 2.; - } - } else { - for (int j = 0; j < sols[i].size(); j++){ - temp[j] = (sols[i][j] + inf_bound.ub()) / 2.; - } - } - if( (_f->eval_vector(temp)[0].ub() < 0.) ){ - zero *current = new zero(sols[i][0], -1); - res->push_back(current); - } else { - zero *current = new zero(sols[i][0], 1); - res->push_back(current); - } - } - return res; -} - -unsigned int integrale::construct_zeros_inf() -{ - _zeros_inf = construct_zeros_intern(_inf_bound); - _zeros_inf->sort(); - return _zeros_inf->size(); -} - -unsigned int integrale::construct_zeros_sup() -{ - _zeros_sup = construct_zeros_intern(_sup_bound); - return _zeros_sup->size(); -} - -ibex::Interval integrale::compute_integral(const double inf, const double sup) -{ - ibex::Interval res(0.0); - double current = inf; - ibex::IntervalVector current_vec(1); - while (current < sup){ - current_vec[0] = ibex::Interval(current, current + __EPS__); - res += _f->eval_vector(current_vec)[0] *__EPS__; - current += __EPS__; - } - - return res; -} - -unsigned int integrale::construct_zeros() -{ - return construct_zeros_inf() + construct_zeros_sup(); -} - -integrale_sol integrale::compute_integral(zero *inf, zero *sup) -{ - ibex::Interval res_inter = compute_integral(inf->get_x().lb(), sup->get_x().ub()); - integrale_sol res(inf, sup, res_inter); - return res; -} - - -std::pair integrale::compute_integral_inf(const ibex::Interval& range) -{ - std::list res_tmp; - if (!_zeros_inf){ - construct_zeros_inf(); - } - if (_zeros_inf->size()){ - std::list::iterator it = _zeros_inf->begin(); - zero *bound_inf = new zero(ibex::Interval(range.lb()), 0); - integrale_sol integ0 = compute_integral(bound_inf, *(it)); - res_tmp.push_back(integ0); - - while (it != _zeros_inf->end()){ - if (*it != _zeros_inf->back()){ - std::list::iterator it_first = it; - it++; - std::list::iterator it_second = it; - integrale_sol integ = compute_integral(*it_first, *it_second); - std::list::iterator it_tmp = res_tmp.begin(); - while (it_tmp != res_tmp.end()){ - *it_tmp = *it_tmp + integ; - it_tmp++; - } - res_tmp.push_back(integ); - - } else { - zero *bound = new zero(ibex::Interval(range.ub()), 0); - integrale_sol integ = compute_integral(*it, bound); - std::list::iterator it_tmp = res_tmp.begin(); - while (it_tmp != res_tmp.end()){ - (*it_tmp) = (*it_tmp) + integ; - it_tmp++; - } - it++; - res_tmp.push_back(integ); - } - } - } else { - zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); - zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); - integrale_sol integ = compute_integral(bound_inf, bound_sup); - res_tmp.push_back(integ); - ibex::Interval nul_inter(0.); - integrale_sol nul(bound_sup, bound_sup, nul_inter); - res_tmp.push_back(nul); - } - if (res_tmp.size() == 0){ - cout << "Erreur, il devrait y avoir au moins une solution" << endl; - exit(EXIT_FAILURE); - } - integrale_sol min = *(res_tmp.begin()); - integrale_sol max = *(res_tmp.begin()); - for (auto elem : res_tmp){ - if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ - min = elem; - } - if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ - max = elem; - } - } - for (auto elem : res_tmp){ - if (elem != min || elem != max){ - if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) - delete elem.get_inf_bound(); - if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) - delete elem.get_sup_bound(); - } - } - return pair(min, max); -} - - -std::pair integrale::compute_integral_sup(const ibex::Interval& range) -{ - std::list res_temp; - if (!_zeros_sup){ - construct_zeros_sup(); - } - if (_zeros_sup->size() > 0){ - std::list::iterator it = _zeros_sup->begin(); - zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); - integrale_sol integ0 = compute_integral((*it), bound_sup); - res_temp.push_back(integ0); - - while (it != _zeros_sup->end()){ - if (*it != _zeros_sup->back()){ - std::list::iterator it_second = it; - it++; - std::list::iterator it_first = it; - integrale_sol integ = compute_integral(*it_first, *it_second); - std::list::iterator it_tmp = res_temp.begin(); - while (it_tmp != res_temp.end()){ - *it_tmp = *it_tmp + integ; - it_tmp++; - } - res_temp.push_back(integ); - - } else { - zero * bound = new zero(ibex::Interval(range.lb()), 0); - integrale_sol integ = compute_integral(bound, *it); - std::list::iterator it_tmp = res_temp.begin(); - while (it_tmp != res_temp.end()){ - *it_tmp = *it_tmp + integ; - it_tmp++; - } - it++; - res_temp.push_back(integ); - } - } - } else { - zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); - zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); - integrale_sol integ = compute_integral(bound_inf, bound_sup); - res_temp.push_back(integ); - ibex::Interval nul_inter(0.); - integrale_sol nul(bound_inf, bound_inf, nul_inter); - res_temp.push_back(nul); - } - if (res_temp.size() == 0){ - cout << "Erreur, il devrait y avoir au moins une solution" << endl; - exit(EXIT_FAILURE); - } - integrale_sol min = *(res_temp.begin()); - integrale_sol max = *(res_temp.begin()); - for (auto elem : res_temp){ - if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ - min = elem; - } - if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ - max = elem; - } - } - for (auto elem : res_temp){ - if (elem != min || elem != max){ - if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) - delete elem.get_inf_bound(); - if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) - delete elem.get_sup_bound(); - } - } - return pair(min, max); -} - -std::pair integrale::compute_integral() -{ - std::pair res = compute_integral_inf(_inf_bound); - zero * inf = res.first.get_sup_bound(); - std::pair res_sup = compute_integral_sup(_sup_bound); - - zero * sup = res_sup.first.get_inf_bound(); - ibex::Interval mid = compute_integral(_inf_bound.ub(), _sup_bound.lb()); - integrale_sol integ_mid(inf, sup, mid); - - res.first+= integ_mid; - res.second+= integ_mid; - res.first+= res_sup.first; - res.second+= res_sup.second; - if (_result == nullptr){ - _result = new std::pair(res); - } else { - *_result = res; - } - return res; -} - - - - -//======================================================================== -std::ostream& operator<<(std::ostream& os, const integrale& inte) -{ - os << "Integral from " << inte.get_inf_bound() << " to " << inte.get_sup_bound() << " of " - << *inte.get_f(); - return os; -} - diff --git a/src/integrale.h b/src/integrale.h deleted file mode 100644 index 196389742d9a9a3408b3e46c0be6746e3bfd1e50..0000000000000000000000000000000000000000 --- a/src/integrale.h +++ /dev/null @@ -1,70 +0,0 @@ -#ifndef _INTEGRALE_ -#define _INTEGRALE_ -#include -#include -#include "zero.h" -#include "integrale_sol.h" - -// #define __DEBUG__ - - - -#define __PREC__ (1e-5) -#define __EPS__ (1e-1) -#define __SHOW_FCT_NAME__ {std::cout << __func__ << std::endl;} - - - -class integrale{ -private: - - ibex::Interval _inf_bound; - ibex::Interval _sup_bound; - - ibex::Function * _f; -// ibex::Variable * _x_var; - - std::list * _zeros_inf; - std::list * _zeros_sup; - - std::pair * _result; - - std::list * construct_zeros_intern(const ibex::Interval & inf_bound); - - integrale_sol compute_integral(zero* inf, zero* sup); - - std::list compute_integral(const double & inf_bound, const ibex::Interval & sup_bound); - - ibex::Interval compute_integral(const double inf_bound, const double sup_bound); - - std::pair compute_integral_inf(const ibex::Interval & range); - - std::pair compute_integral_sup(const ibex::Interval & range); - -public: - //CONSTRUCTORS - integrale(const ibex::Interval& inf, const ibex::Interval& sup, ibex::Function * f); - ~integrale(); - - //GETTERS - inline ibex::Interval get_inf_bound() const {return _inf_bound;}; - inline ibex::Interval get_sup_bound() const {return _sup_bound;}; - inline ibex::Function * get_f() const {return _f;}; -// inline ibex::Variable * get_x_var() const {return _x_var;}; - inline std::list * get_zeros_inf(){return _zeros_inf;}; - inline std::list * get_zeros_sup(){return _zeros_sup;}; - inline std::pair * get_result(){return _result;}; - - unsigned int construct_zeros_inf(); - - unsigned int construct_zeros_sup(); - - unsigned int construct_zeros(); - - std::pair compute_integral(); -}; - - -std::ostream& operator<<(std::ostream& os, const integrale& inte); - -#endif //_INTEGRALE_ diff --git a/src/integrale_sol.cpp b/src/integrale_sol.cpp deleted file mode 100644 index f4f20fc8eb029cefe58debe57bdd62ca00dffc23..0000000000000000000000000000000000000000 --- a/src/integrale_sol.cpp +++ /dev/null @@ -1,76 +0,0 @@ -/* - * - * Copyright (C) 2019 - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include "integrale_sol.h" - -integrale_sol::integrale_sol(zero* inf, zero* sup, const ibex::Interval& res) : _res_integrale(res) -{ - _inf = (inf); - _sup = (sup); - _min_candidate = false; - _max_candidate = false; - if (res.ub() < 0.){ - _min_candidate = true; - } - if (res.lb() > 0.){ - _max_candidate = true; - } -} - -integrale_sol::~integrale_sol() -{ -// if (_inf) -// delete _inf; -// if (_sup) -// delete _sup; -} - - -integrale_sol & integrale_sol::operator+=(const integrale_sol &rhs){ - *this = *this + rhs; - return *this; -} - -integrale_sol operator+(const integrale_sol& lhs, const integrale_sol& rhs) -{ - - if (lhs.get_sup_bound() == rhs.get_inf_bound()){ - return integrale_sol(lhs.get_inf_bound(), rhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); - } else { - return integrale_sol(rhs.get_inf_bound(), lhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); - } - -} - -std::ostream & operator<<(std::ostream& os, const integrale_sol& inte) -{ - os << "Integral from (" << *(inte.get_inf_bound()) << ") to (" << *inte.get_sup_bound() << ") : " << inte.get_res_integrale(); - return os; -} - -bool operator==(const integrale_sol& lhs, const integrale_sol& rhs) -{ - return (lhs.get_inf_bound() == rhs.get_inf_bound() && lhs.get_sup_bound() == rhs.get_sup_bound()); -} - -bool operator!=(const integrale_sol& lhs, const integrale_sol& rhs) -{ - return !(lhs == rhs); -} - - diff --git a/src/integrale_sol.h b/src/integrale_sol.h deleted file mode 100644 index c30247ee00403467bedc1938dd1542f0942f0f31..0000000000000000000000000000000000000000 --- a/src/integrale_sol.h +++ /dev/null @@ -1,57 +0,0 @@ -/* - * - * Copyright (C) 2019 - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#ifndef INTEGRALE_SOL_H -#define INTEGRALE_SOL_H - -#include "zero.h" - -/** - * @todo write docs - */ -class integrale_sol -{ -private: - zero * _inf; - zero * _sup; - ibex::Interval _res_integrale; - bool _min_candidate; - bool _max_candidate; - -public: - integrale_sol(zero *inf, zero *sup, const ibex::Interval & res); - ~integrale_sol(); - - inline zero * get_inf_bound() const {return _inf;}; - inline zero * get_sup_bound() const {return _sup;}; - inline ibex::Interval get_res_integrale() const {return _res_integrale;}; - inline bool is_min_candidate(){return _min_candidate;}; - inline bool is_max_candidate(){return _max_candidate;}; - void set_min_candidate(bool min){_min_candidate = min;}; - void set_max_candidate(bool max){_max_candidate = max;}; - - - integrale_sol & operator+=(const integrale_sol &rhs); -}; - -integrale_sol operator+(const integrale_sol &, const integrale_sol &); -bool operator==(const integrale_sol &, const integrale_sol &); -bool operator!=(const integrale_sol &, const integrale_sol &); -std::ostream& operator<<(std::ostream& os, const integrale_sol& inte_sol); - -#endif // INTEGRALE_SOL_H diff --git a/src/integralii.cpp b/src/integralii.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b09222ca32902b42c0c9c234b3939e723de35290 --- /dev/null +++ b/src/integralii.cpp @@ -0,0 +1,176 @@ +/* + * + * Copyright (C) 2019 mullier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "integralii.h" + + +std::list& IntegralII::compute_list_of_zeros(ibex::Interval searchdomain) +{ + //ATTENTION pas de preuve d'unicité + std::list* res = new std::list(); + + ibex::CtcNewton newton(*_f); + + ibex::CtcFwdBwd fb (*_f); + ibex::RoundRobin bisector(0); + + ibex::CellStack buff; + + ibex::CtcCompo comp( fb, newton); + + ibex::CtcFixPoint fix(comp, 1e-7); + + ibex::Solver s(fix, bisector, buff); + + + ibex::IntervalVector inf(1); + inf[0] = searchdomain; + std::vector sols = s.solve(inf); + + unsigned int i; + for (i = 0; i < sols.size(); i++){ + //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 + ibex::Interval * current = new ibex::Interval(sols[i][0]); + res->push_back(current); + } + return *res; +} + +std::list& IntegralII::compute_list_of_zeros_inf() +{ + std::list& res = compute_list_of_zeros(_inf_endpoint); + res.sort(); + return res; +} + +std::list& IntegralII::compute_list_of_zeros_sup() +{ + return compute_list_of_zeros(_sup_endpoint); +} + + +void IntegralII::compute_integral() +{ + + std::list res_temp; + std::list::iterator it_sols; + ibex::Interval *current = new ibex::Interval(_inf_endpoint.lb()); + IntegralRR min, max; + //======================================================================= + //Computation of the integrals in the _inf_endpoint + std::list list_of_zeros_inf = compute_list_of_zeros_inf(); + if (list_of_zeros_inf.size()>0){ + std::list::iterator it_zeros_inf = list_of_zeros_inf.begin(); + while (it_zeros_inf != list_of_zeros_inf.end()){//calcul des intégrales entre chaque zéro de f + IntegralRR current_integral(current->mid(), (*it_zeros_inf)->mid(), _f); + it_sols = res_temp.begin(); + while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales + *it_sols += current_integral; + it_sols++; + } + res_temp.push_back(current_integral); + current = *it_zeros_inf; + it_zeros_inf++; + } + //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle + IntegralRR current_integral(current->mid(), _inf_endpoint.ub(), _f); + it_sols = res_temp.begin(); + while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales + *it_sols += current_integral; + it_sols++; + } + res_temp.push_back(current_integral); + } else { + IntegralRR whole_integral(_inf_endpoint.lb(), _inf_endpoint.ub(), _f); + res_temp.push_back(whole_integral); + res_temp.push_back(IntegralRR(_inf_endpoint.ub(), _inf_endpoint.ub(), _f, 0.)); + } + it_sols = res_temp.begin(); + min = *it_sols; + max = *it_sols; + it_sols++; + while (it_sols != res_temp.end()){//On recherche le min et le max + if (max < *it_sols){ + max = *it_sols; + } + if (*it_sols < min){ + min = *it_sols; + } + it_sols ++; + } + delete current; + + //Calcul de l'intégrale entre les intervalles + IntegralRR integrale_milieu(_inf_endpoint.ub(), _sup_endpoint.lb(), _f); + _solution.first = min + integrale_milieu; + _solution.second = max + integrale_milieu; + + //======================================================================= + current = new ibex::Interval(_sup_endpoint.ub()); + res_temp.clear(); + std::list list_of_zeros_sup = compute_list_of_zeros_sup(); + if (list_of_zeros_sup.size()>0){ + std::list::iterator it_zeros_sup = list_of_zeros_sup.begin(); + while (it_zeros_sup != list_of_zeros_sup.end()){//calcul des intégrales entre chaque zéro de f + IntegralRR current_integral((*it_zeros_sup)->mid(), current->mid(), _f); + it_sols = res_temp.begin(); + while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales + *it_sols += current_integral; + it_sols++; + } + res_temp.push_back(current_integral); + current = *it_zeros_sup; + it_zeros_sup++; + } + //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle + IntegralRR current_integral(_sup_endpoint.lb(), current->mid(), _f); + it_sols = res_temp.begin(); + while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales + *it_sols += current_integral; + it_sols++; + } + res_temp.push_back(current_integral); + } else { + IntegralRR whole_integral(_sup_endpoint.lb(), _sup_endpoint.ub(), _f); + res_temp.push_back(whole_integral); + res_temp.push_back(IntegralRR(_sup_endpoint.lb(), _sup_endpoint.lb(), _f, 0.)); + } + it_sols = res_temp.begin(); + min = *it_sols; + max = *it_sols; + it_sols++; + while (it_sols != res_temp.end()){//On recherche le min et le max + if (max < *it_sols){ + max = *it_sols; + } + if (*it_sols < min){ + min = *it_sols; + } + it_sols ++; + } + + _solution.first += min; + _solution.second += max; + +} + +std::ostream& operator<<(std::ostream& os, const IntegralII& integral){ + os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); + os << " of the function " << *integral.get_f() << " : " << endl << "minimum : " << integral.get_solution().first << endl << "maximum : " << integral.get_solution().second; + return os; +} diff --git a/src/integralii.h b/src/integralii.h new file mode 100644 index 0000000000000000000000000000000000000000..72ae2ee70911f2a2c1883f598f961ef63a1634fa --- /dev/null +++ b/src/integralii.h @@ -0,0 +1,60 @@ +/* + * + * Copyright (C) 2019 mullier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef INTEGRALII_H +#define INTEGRALII_H + +#include +#include "integralrr.h" + +/** + * @todo write docs + */ +class IntegralII +{ +private: + ibex::Interval _inf_endpoint; + ibex::Interval _sup_endpoint; + + ibex::Function * _f; + + std::pair _solution; + + std::list& compute_list_of_zeros(ibex::Interval searchdomain); + +public: + IntegralII(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function * f) : _inf_endpoint(inf), _sup_endpoint(sup), _f(f){}; + + ibex::Interval get_inf_endpoint() const {return _inf_endpoint;}; + ibex::Interval get_sup_endpoint() const {return _sup_endpoint;}; + + ibex::Function * get_f() const {return _f;}; + + std::pair get_solution() const {return _solution;}; + + std::list& compute_list_of_zeros_inf(); + + std::list& compute_list_of_zeros_sup(); + + void compute_integral(); + +}; + +std::ostream& operator<<(std::ostream& os, const IntegralII& integral); + +#endif // INTEGRALEII_H diff --git a/src/integralrr.cpp b/src/integralrr.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8f628f8c09d20328b199f0ffa5c80cd6d4cab5b7 --- /dev/null +++ b/src/integralrr.cpp @@ -0,0 +1,69 @@ +/* + * + * Copyright (C) 2019 mullier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "integralrr.h" + +IntegralRR::IntegralRR(double inf_endpoint, double sup_endpoint, ibex::Function* f) : _inf_endpoint(inf_endpoint), _sup_endpoint(sup_endpoint), _f(f), _solution(0.) +{ + compute_integral(); +} + +void IntegralRR::compute_integral() +{ + ibex::Interval res(0.0); + double current = _inf_endpoint; + ibex::IntervalVector current_vec(1); + while (current < _sup_endpoint){ + current_vec[0] = ibex::Interval(current, current + __EPS__); + res += _f->eval_vector(current_vec)[0] *__EPS__; + current += __EPS__; + } + current_vec[0] = ibex::Interval(current - __EPS__, _sup_endpoint); + res += _f->eval_vector(current_vec)[0] * (_sup_endpoint - current + __EPS__); + _solution = res; +} + + +IntegralRR& IntegralRR::operator+=(const IntegralRR& other) +{ + *this = *this + other; + return *this; +} + + +IntegralRR operator+(const IntegralRR& lhs, const IntegralRR& rhs) +{ + if (lhs.get_sup_endpoint() == rhs.get_inf_endpoint()){ + return IntegralRR(lhs.get_inf_endpoint(), rhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); + } else { + return IntegralRR(rhs.get_inf_endpoint(), lhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); + } + +} + +std::ostream & operator<<(std::ostream& os, const IntegralRR& integral) +{ + os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); + os << " of the function " << *integral.get_f() << " : " << integral.get_solution(); + return os; +} + +bool operator<(const IntegralRR& lhs, const IntegralRR& rhs) +{ + return lhs.get_solution().ub() < rhs.get_solution().lb(); +} diff --git a/src/integralrr.h b/src/integralrr.h new file mode 100644 index 0000000000000000000000000000000000000000..8bcd766195d9ad21c260b4b082cec3e52a25d6ff --- /dev/null +++ b/src/integralrr.h @@ -0,0 +1,67 @@ +/* + * + * Copyright (C) 2019 mullier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef INTEGRALRR_H +#define INTEGRALRR_H + +#include + +#define __EPS__ (1e-4) + +/** + * @todo write docs + */ +class IntegralRR +{ +private: + double _inf_endpoint; + double _sup_endpoint; + + ibex::Function * _f; + + ibex::Interval _solution; + +public: + //CONSTRUCTORS + IntegralRR(double, double, ibex::Function *); + IntegralRR(double inf, double sup, ibex::Function * f, const ibex::Interval & sol) : + _inf_endpoint(inf), _sup_endpoint(sup), _f(f), _solution(sol){}; + IntegralRR(){_f = nullptr;}; + IntegralRR(const IntegralRR & other) : IntegralRR(other.get_inf_endpoint(), other.get_sup_endpoint(), other.get_f(), other.get_solution()){}; + + //GETTERS + inline double get_inf_endpoint() const {return _inf_endpoint;}; + inline double get_sup_endpoint() const {return _sup_endpoint;}; + inline ibex::Function * get_f() const {return _f;}; + inline ibex::Interval get_solution() const {return _solution;}; + + //SETTERS + void set_solution(const ibex::Interval & res){_solution = res;}; + + void compute_integral(); + + IntegralRR& operator+=(const IntegralRR& other); +}; + +IntegralRR operator+(const IntegralRR& lhs, const IntegralRR& rhs); + +bool operator<(const IntegralRR &, const IntegralRR &); + +std::ostream& operator<<(std::ostream& os, const IntegralRR& integral); + +#endif // INTEGRALRR_H diff --git a/src/main.cpp b/src/main.cpp index 70498dd3303c2533dcfeed69f55669f8a86e61e7..702301d01ce65f33c1f043b71b53097875d6a8a6 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,33 +1,17 @@ -#include -#include "integrale.h" +#include -using namespace std; +#include "integralii.h" - - -void example(){ +int main(){ cout << setprecision(15); ibex::Variable x(1); - ibex::Function f(x, -15. + x+5.*ibex::sin(x)); - ibex::Interval inf(0., 5.5); - ibex::Interval sup(10, 15.); - - integrale example(inf, sup, &f); - - cout << example << endl; - example.construct_zeros(); - for (auto elem : *example.get_zeros_sup()){ - cout << elem << endl; - cout << *elem << endl << endl; - } - example.compute_integral(); - cout << "minimum: " << example.get_result()->first << endl; - cout << "maximum: " << example.get_result()->second << endl; + ibex::Function f(x, ibex::sin(x)); + ibex::Interval inf_inter(0, 5); + ibex::Interval sup_inter(9, 10); -} - - -int main(/*int argc, char * argv[]*/){ - example(); + IntegralII integ_inter(inf_inter, sup_inter, &f); + integ_inter.compute_integral(); + cout << "Result: " << endl; + cout << integ_inter << endl; return EXIT_SUCCESS; } diff --git a/src/zero.cpp b/src/zero.cpp deleted file mode 100644 index 64119f83bf29b7cb3534ecd36f637921d0a63dfd..0000000000000000000000000000000000000000 --- a/src/zero.cpp +++ /dev/null @@ -1,33 +0,0 @@ -/* - * - * Copyright (C) 2019 - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include "zero.h" - -ostream& operator<<(ostream& os, const zero& z) -{ - if (z.is_non_zero()) - os << "boundary: " << z.get_x(); - else { - os << "Zero of the function: " << z.get_x() << ", derivative is "; - if (z.is_neg()) - os << " negative"; - else - os << " positive"; - } - return os; -} diff --git a/src/zero.h b/src/zero.h deleted file mode 100644 index b4fe88b17a442abda35600b27b07e25fda87ce2f..0000000000000000000000000000000000000000 --- a/src/zero.h +++ /dev/null @@ -1,55 +0,0 @@ -/* - * - * Copyright (C) 2019 - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#ifndef ZERO_H -#define ZERO_H - -#include - -/** - * @todo write docs - */ -class zero -{ -private: - ibex::Interval _x; - int _slope; - -public: - zero(const zero &z) : _x(z.get_x()), _slope(z.get_slope()){}; - zero(ibex::Interval x, int slope) : _x(x), _slope(slope) {}; - - inline ibex::Interval get_x() const {return _x;}; - inline int get_slope() const {return _slope;}; - - inline bool is_neg() const {return _slope < 0;}; - inline bool is_pos() const {return _slope > 0;}; - inline bool is_non_zero() const {return _slope == 0;}; -}; - -inline bool operator< (const zero& lhs, const zero& rhs){ return lhs.get_x().mid() < rhs.get_x().mid(); } -inline bool operator> (const zero& lhs, const zero& rhs){ return rhs < lhs; } -inline bool operator<=(const zero& lhs, const zero& rhs){ return !(lhs > rhs); } -inline bool operator>=(const zero& lhs, const zero& rhs){ return !(lhs < rhs); } - -inline bool operator==(const zero& lhs, const zero& rhs){ return lhs.get_x() == rhs.get_x(); } -inline bool operator!=(const zero& lhs, const zero& rhs){ return !(lhs == rhs); } - -ostream& operator<<(ostream& os, const zero& z); - -#endif // ZERO_H