# TP3: Méthodes de points fixes
## Exercice 1
### Faire des graphs pour représenter les fonctions et les points fixes

```{r}
g1 <- function(x) {
    (x^2 - 1) / 3
}
g1p <- function(x){abs(2*x/3)}
curve(g1, from = -1, to = 4)
abline(0, 1)
points((3-sqrt(13))/2, (3-sqrt(13))/2, col="red")
points((3+sqrt(13))/2, (3+sqrt(13))/2, col="red")
```

## Exercice 2
### Question 2

```{r}
pfixe <- function(x0, eps, fun, Nmax){
    ## Initialisation
    n <- 0
    xn <- x0
    xnp1 <- fun(xn)
    ## Iterations, on s'arrete quand l'erreur est atteinte
    while(abs(xn-xnp1) > eps & n < Nmax){
        xn <- xnp1
        xnp1 <- fun(xn)
        n <- n + 1
    }
    ## On renvoit la valeur de xn, et le nombre d'itérations
    return(list("x"=xnp1,"iter"=n+1))
}
```

### Question 3
#### Remarque, nous avons une erreur pour la fonction g3 parce que l'algorithme diverge et la valeur est trop grande.
```{r}
f <- function(x) {x^2-2*x-3}
g1 <- function(x){sqrt(2*x+3)} # On va converger
g2 <- function(x){3/(x-2)} # On va converger
g3 <- function(x){(x^2-3)/2} # On va diverger

x0 <- 4
eps <- 1e-8
Nmax <- 100

resg1 <- pfixe(x0=x0,
            eps=eps,
            fun=g1,
            Nmax=Nmax)

resg2 <- pfixe(x0=x0,
            eps=eps,
            fun=g2,
            Nmax=Nmax)

resg3 <- pfixe(x0=x0,
        eps=eps,
        fun=g3,
        Nmax=Nmax)

resg2
```

### Question 4

```{r}
pfixe2 <- function(x0, eps, fun, Nmax){
    ## Initialisation
    n <- 0
    xn <- x0
    xnp1 <- fun(xn)
    ## Initialisation d'une liste pour garder les valeurs successives de xn a chaque intération
    xnp <- c()
    while (abs(xn-xnp1) > eps & n < Nmax){
        xn <- xnp1
        xnp1 <- fun(xn)
        n <- n + 1
        xnp[n] <- xnp1
    }
    return(list("x"=c(fun(x0),xnp),
    "n"=n+1))
}
resg1_2 <- pfixe2(x0=x0,
            eps=eps,
            fun=g1,
            Nmax=Nmax)
resg2_2 <- pfixe2(x0=x0,
            eps=eps,
            fun=g2,
            Nmax=Nmax)
resg3_2 <- pfixe2(x0=x0,
            eps=eps,
            fun=g3,
            Nmax=2)

```

#### Remarque, on utilise la formule théorique qu'on a calculé dans l'exercice 1.a

```{r}
iterpfixe <- function(x0, x1, k, eps) {
    ceiling(log(eps*(1-k)/abs(x1-x0))/log(k))
}
k <- 0.4
x1 <- g2(x0)
iterpfixe(x0, x1, k, eps)
```

#### Remarque, dans les exemples suivants, on plot successivement la courbe dont les points sont (xn, xn+1)

```{r}
g1x <- rep(resg1_2$x, each=2)
xg1 <- c(x0, head(g1x,-1))
plot(xg1,g1x,type="s")
curve(g1,add=TRUE, col="blue")
abline(0,1,col="red")
```

```{r}
g2x <- rep(resg2_2$x, each=2)
xg2 <- c(x0, head(g2x,-1))
plot(xg2,g2x,type="s")
curve(g2,add=TRUE, col="blue")
abline(0,1,col="red")
```

```{r}
g3x <- rep(resg3_2$x, each=2)
xg3 <- c(x0, head(g3x,-1))
plot(xg3,g3x,type="s")
curve(g3,add=TRUE, col="blue")
abline(0,1,col="red")
```
```{r}
res3 <- pfixe2(x0=-0.5, eps=eps, fun=g3, Nmax=200)
g3x <- rep(res3$x, each=2)
xg3 <- c(x0, head(g3x,-1))
plot(xg3,g3x,type="s")
curve(g3, add=T, col="blue")
abline(0,1,col="red")
```

```{r}
plot(res3$x,type="l")
```

### Question 5

```{r}
g5 <- function(x){sqrt(1+x)}
x0 <- 1.5
resg5 <- pfixe2(x0=x0, eps=eps, fun=g5, Nmax=Nmax)

en <- abs(resg5$x-resg5$x[length(resg5$x)])

plot(en, type="l")
```

```{r}
x_en <- abs(en[1:(length(en)-2)])
y_en <- abs(en[2:(length(en)-1)])
plot(x_en, y_en, log="xy", type="l", col = "green")
C = median(abs(en[2:length(en)]/en[1:(length(en)-1)]))
p = 1
lines(x_en, p*x_en, col="red", type="l")
C
```

### Question 6

```{r}

gN <- function(x){x-(x^2-2)/(2*x)}
resN <- pfixe2(x0=x0,
eps=eps,
fun=gN,
Nmax=Nmax)
enN <- resN$x-sqrt(2)
xN <- log(abs(enN[1:(length(enN)-2)]))
yN <- log(abs(enN[2:(length(enN)-1)]))
plot(xN, yN,
type="l")
p = 2
lines(xN,
p*xN,col="red",type="l")
C = median(abs(enN[2:length(enN)]/(enN[1:(length(enN)-1)])^2))
C
```

## Exercice 4
### Question 1

```{r}

Newton2d <- function(x0, eps=.Machine$double.eps^(1/2), fun, J, Nmax){
    ## Vecnorm nous aide à calculer notre erreur en grande dimension
    vecnorm <- function(x) { sqrt(sum(x^2)) }
    ## Initialisation
    n <- 0
    err <- eps + 1
    x <- x0
    while (err > eps && n < Nmax){
        n <- n + 1
        Fx <- fun(x)
        Jx <- J(x)
        ## On résoud le system Jx*z = Fx
        z <- solve(Jx, Fx)
        x <- x - z
        err <- vecnorm(z)
    }
    Fx <- vecnorm(fun(x))
    if (n > Nmax && err > eps) {
        cat("Erreur trop grande apres le nombre d'itération maximum, absence de convergence")
    }
    return(list(zero=x, fnorm=Fx, niter=n))
}

```

### Question 2

```{r}
## Systeme tel que f(x) = 0
F1 <- function(x) c(3*x[1]^2 - x[2]^2, 3*x[1]*x[2]-x[1]^3-1)
## Calcul de la jacobienne
JF1 <- function(x) matrix( c(6*x[1], -2*x[2], 3*x[2]-3*x[1]^2, 3*x[1]), 2, 2, byrow = T)
## x0
X0 <- c(1,1)
eps <- 1e-8
Nmax <- 100
Newton2d(x0 = X0, eps = eps, fun = F1, J = JF1, Nmax = Nmax)
```