integralrr.cpp 2.38 KB
Newer Older
Olivier Mullier's avatar
Olivier Mullier committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
/*
 * <one line to give the library's name and an idea of what it does.>
 * Copyright (C) 2019  mullier <email>
 *
 * 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 <http://www.gnu.org/licenses/>.
 */

#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();
}