6663b6c9
 adorian
projet complet av...
|
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
#include <poincare/derivative.h>
#include <poincare/symbol.h>
#include <poincare/simplification_engine.h>
#include <poincare/undefined.h>
#include <cmath>
extern "C" {
#include <assert.h>
#include <float.h>
}
namespace Poincare {
Expression::Type Derivative::type() const {
return Type::Derivative;
}
Expression * Derivative::clone() const {
Derivative * a = new Derivative(m_operands, true);
return a;
}
int Derivative::polynomialDegree(char symbolName) const {
if (symbolName == 'x') {
if (operand(1)->polynomialDegree(symbolName) != 0) {
return -1;
}
return 0;
}
return Expression::polynomialDegree(symbolName);
}
Expression * Derivative::shallowReduce(Context& context, AngleUnit angleUnit) {
Expression * e = Expression::shallowReduce(context, angleUnit);
if (e != this) {
return e;
}
#if MATRIX_EXACT_REDUCING
if (operand(0)->type() == Type::Matrix || operand(1)->type() == Type::Matrix) {
return replaceWith(new Undefined(), true);
}
#endif
// TODO: to be implemented diff(+) -> +diff() etc
return this;
}
template<typename T>
Complex<T> * Derivative::templatedApproximate(Context& context, AngleUnit angleUnit) const {
static T min = sizeof(T) == sizeof(double) ? DBL_MIN : FLT_MIN;
static T epsilon = sizeof(T) == sizeof(double) ? DBL_EPSILON : FLT_EPSILON;
Evaluation<T> * xInput = operand(1)->privateApproximate(T(), context, angleUnit);
T x = xInput->toScalar();
delete xInput;
T functionValue = operand(0)->approximateWithValueForSymbol('x', x, context, angleUnit);
// No complex/matrix version of Derivative
if (std::isnan(x) || std::isnan(functionValue)) {
return new Complex<T>(Complex<T>::Undefined());
}
T error, result;
T h = k_minInitialRate;
do {
result = riddersApproximation(context, angleUnit, x, h, &error);
h /= 10.0;
} while ((std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) && h >= epsilon);
/* if the error is too big regarding the value, do not return the answer */
if (std::fabs(error/result) > k_maxErrorRateOnApproximation || std::isnan(error)) {
return new Complex<T>(Complex<T>::Undefined());
}
if (std::fabs(error) < min) {
return new Complex<T>(result);
}
error = std::pow((T)10, std::floor(std::log10(std::fabs(error)))+2);
return new Complex<T>(std::round(result/error)*error);
}
template<typename T>
T Derivative::growthRateAroundAbscissa(T x, T h, Context & context, AngleUnit angleUnit) const {
T expressionPlus = operand(0)->approximateWithValueForSymbol('x', x+h, context, angleUnit);
T expressionMinus = operand(0)->approximateWithValueForSymbol('x', x-h, context, angleUnit);
return (expressionPlus - expressionMinus)/(2*h);
}
template<typename T>
T Derivative::riddersApproximation(Context & context, AngleUnit angleUnit, T x, T h, T * error) const {
/* Ridders' Algorithm
* Blibliography:
* - Ridders, C.J.F. 1982, Advances in Engineering Software, vol. 4, no. 2,
* pp. 75–76. */
*error = sizeof(T) == sizeof(float) ? FLT_MAX : DBL_MAX;
// Initialize hh
assert(h != 0.0);
/* Make hh an exactly representable number */
volatile T temp = x+h;
T hh = temp - x;
/* a is matrix storing the function extrapolations for different stepsizes at
* different order */
T a[10][10];
for (int i = 0; i < 10; i++) {
for (int j = 0; j < 10; j++) {
a[i][j] = 1;
}
}
a[0][0] = growthRateAroundAbscissa(x, hh, context, angleUnit);
T ans = 0;
T errt = 0;
/* Loop on i: change the step size */
for (int i = 1; i < 10; i++) {
hh /= k_rateStepSize;
/* Make hh an exactly representable number */
volatile T temp = x+hh;
hh = temp - x;
a[0][i] = growthRateAroundAbscissa(x, hh, context, angleUnit);
T fac = k_rateStepSize*k_rateStepSize;
/* Loop on j: compute extrapolation for several orders */
for (int j = 1; j < 10; j++) {
a[j][i] = (a[j-1][i]*fac-a[j-1][i-1])/(fac-1);
fac = k_rateStepSize*k_rateStepSize*fac;
errt = std::fabs(a[j][i]-a[j-1][i]) > std::fabs(a[j][i]-a[j-1][i-1]) ? std::fabs(a[j][i]-a[j-1][i]) : std::fabs(a[j][i]-a[j-1][i-1]);
/* Update error and answer if error decreases */
if (errt < *error) {
*error = errt;
ans = a[j][i];
}
}
/* If higher extrapolation order significantly increases the error, return
* early */
if (std::fabs(a[i][i]-a[i-1][i-1]) > 2*(*error)) {
break;
}
}
return ans;
}
}
|