-
Notifications
You must be signed in to change notification settings - Fork 0
/
brent.c
110 lines (99 loc) · 2.17 KB
/
brent.c
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
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "brent.h"
#define numerador(x0,fx1,fm) x0*fx1*fxm
#define denominador(fx0,fx1,fxm) (fx0 - fx1 )*( fx0 - fxm )
#define parcela(x0,fx0,fx1,fxm) numerador(x0,fx1,fxm)/denominador(fx0,fx1,fxm)
#define calc_inv(fx0,fx1,fxm,x0,x1,xm) (parcela(x0,fx0,fx1,fxm) + parcela(x1,fx1,fx0,fxm) +parcela(xm,fxm,fx0,fxm))
#define secante(a,b,fa,fb) (b-(fb*((b-a)/(fb-fa))))
double brent (double x0, double x1, int p, double(*f) (double x)){
double precision = 5 * pow((double)10,-p);
double xm =0,fx0,fx1,fxm,x2,fx2,erro,modulo,mod1,trap;
int usou = 0;
xm = (x0 +x1)/2;
fxm = f(xm);
trap = fabs(fxm);
while(fabs(fxm) > precision ){
//printf("__fxm:%g___precision:%g\n",fabs(fxm),precision);
if(usou && (fabs(trap - fxm) < precision) ){
//Não convergirá
printf("Nao convergiu a f(x)=0\n");
break;
}
else{
trap = fxm;
}
puts("passo");
getchar();
fx0 = f(x0);
fx1 = f(x1);
usou = 0;
/////IQI
x2 =(double) calc_inv(fx0,fx1,fxm,x0,x1,xm) ;
erro = fabs(f(xm));
fx2= f(x2);
//printf("__x:%g___f():%g\n",x2,fx2);
//Testando f(IQI)
if (fabs(f(x2)) < erro){
modulo = fabs((x0-x1));
if(fx0 *fx2 < 0){
mod1 = fabs((x0-x2));
if(mod1/modulo <=0.5){
x0 = x2;
usou =1;
puts("A");
}
}
else{
mod1 = fabs(x1-x2);
if(mod1/modulo <=0.5){
x1=x2;
usou = 1;
puts("A");
}
}
}
if (!usou){
x2 = (double) secante(x0,x1,fx0,fx1);
//Testando f(secante)
if (f(x2) < erro){
modulo = fabs((x0-x1));
if(fx0 *fx2 < 0){
mod1 = fabs((x0-x2));
if(mod1/modulo <=0.5){
x1 = x2;
usou =1;
puts("B");
}
}
else{
mod1 = fabs(x1-x2);
if(mod1/modulo <=0.5){
x0=x2;
usou = 1;
puts("B");
}
}
}
}
if(!usou){
//usa xm
if(fx0 *fxm < 0){
x1=xm;
usou = 1;
puts("C");
}
else{
x0=xm;
usou = 1;
puts("C");
}
}
xm = (x0 +x1)/2;
fxm = f(xm);
}
//printf("__fxm:%g___precision:%g\n",fxm,precision);
//printf("__x1:%g___x2:%g\n",x0,x1);
return xm;
}