Pour \(n\geq 1\), (par construction) on a \(b_n-a_n=\frac{b-a}{2^{n-1}}\). Comme \(x_n=(a_n+b_n)/2\), pour tout \(n\geq 1\) et \(x^*\in [a_n,b_n]\), on a \(|x_n-x^*| \leq (b_n-a_n)/2=(b-a)/2^n\). La suite \((x_n)_{n\geq 1}\) converge vers \(x^*\).
La deuxième question consiste à déterminer à partir de quel rang \(n_0\) a-t-on \(|x_n-x^*|\leq \varepsilon\), où \(\varepsilon >0\). i.e \(\forall \varepsilon >0, \exists~n>n_0\) tel que \(|x_n-x_0|<\varepsilon\). Il suffit d’avoir \((b-a)/2^n \leq \varepsilon\). On peut prendre, \(n_0 = E( \log_2((b-a)/\varepsilon) ) + 1.\) (\(E\) partie entière). \(\forall \epsilon > 0\), si \(n\geq n_0\), alors \(|x_n-x^*|\leq \epsilon\).
f <- function(x){
y = exp(x)-4*x
return(y)
#x^3-3*x+1
}
curve(f, from = 0, to = 3)
abline(h=0)
dicho <- function(a, b, eps, f, Nmax){
#
# methode de dicho pour resoudre f(x)=0
#
# input :
# - a borne inf de l'intervale de recherche
# - b borne sup //
# - eps precision sur la solution
# - f fonction pour laquelle on cherche un zero
# - Nmax nombre max d'iteration
# ouput :
# - xn la solution
# - n le nb d'iteration
n <- 0
an <- a
bn <- b
xn <- (an + bn)/2
while (abs(bn-an) > eps & n < Nmax){
if (f(an)*f(xn) < 0){
bn <- xn
}else{
an <- xn
}
xn <- (an + bn)/2
n <- n + 1
}
return(list("x" = xn,
"n" = n))
}
eps <- 1e-8
Nmax <- 100
a <- 1
b <- 2.5
res <- dicho(a=a, b=b, eps=eps, f=f, Nmax=Nmax)
res
## $x
## [1] 2.153292
##
## $n
## [1] 28
curve(f, from = 1, to = 2.5)
abline(h=0)
points(res$x, f(res$x), col="red", pch=19)
iterdicho <- function(a, b, eps){
L <- b - a
y <- floor(log2(L/eps)+ 1)
y
}
print(res$n)
## [1] 28
print(iterdicho(a, b , eps))
## [1] 28
dicho2 <- function(a, b, eps, f, Nmax){
if(a >= b){stop("a must be < b")}
n <- 1
an <- a
bn <- b
xn <- (an + bn)/2
xn_n <- c()
xn_n[1] <- xn
while (abs(bn-an) > eps & n < Nmax){
if (f(an)*f(xn) < 0){
bn <- xn
}else{
an <- xn
}
xn <- (an + bn)/2
xn_n[n+1] <- xn
n <- n + 1
}
if(n == Nmax){
stop("n=Nmax")
}else{
message("la methode converge après ", n-1, " itérations")
}
return(list("x" = xn,
"n" = n,
"xn_n" = xn_n))
}
res2 <- dicho2(a, b, eps, f, Nmax)
## la methode converge après 28 itérations
#xn <- c( (a+b)/2, res2$xn_n)
err <- abs(res2$xn_n-res2$x)
plot(1:res2$n, (b-a)/2^(1:res2$n),
type="l", log = "y", col="red")
lines(1:res2$n, err, type="l")
On pose \(e_n=|x^*-x_{n}|\), on a \(|e_{n+1}|\leq \frac{1}{2}|e_n|\), à chaque itération on divise par deux l’erreur. La méthode est d’ordre 1.
err[length(err)-1]/err[(length(err)-2)]
## [1] 0.3333333
Estimation par régression linéaire
xn <- log(err[1:(length(err)-2)])
yn <- log(err[2:(length(err)-1)])
reg <- lm(yn~xn)
hat_p <- coef(reg)["xn"]
hat_p # estimation de l'ordre
## xn
## 0.9586508
hat_C <- exp(coef(reg)["(Intercept)"])
hat_C # estimation de la vitesse
## (Intercept)
## 0.3255163
plot(xn,yn)
abline(reg)
k0 <- 20
k1 <- 75
alpha <- 1.5
D <- function(p){
y<- k0*(1+1/p)
}
O <- function(p){
y <- k1*p^alpha
}
OmD <- function(p){
y <- k1*p^alpha - k0*(1+1/p)
}
curve(D, from = 0, to = 1)
curve(O, from = 0, to = 1, add = T)
pstar <- dicho(a=0.5, b=1, eps=eps, f=OmD, Nmax=Nmax)
pstar$x
## [1] 0.7346332
pstar$n
## [1] 26