Draw fractals from root finding iteration in R

class: center, middle, inverse, title-slide

# Draw fractals from root finding iteration in R
## LA R users group: April Meeting
### Keren Xu
### PhD student in Epidemiology at USC
### 2020/04/23

---

# Overview   
### 1. Root-finding Algorithm

### 2. Newton’s Method

### 3. Secant Method

### 4. Fractals

### 5. Getting Started - Creating your fractals

---
# Root-finding Algorithm

`\(f(a) = 0\)`

`\(f(x) = g(x)\)` ➡️ `\(h(x) = f(x) – g(x)\)`

.center[
<img src="root_def.png" width="50%" height="50%" />
]

---

background-image: url(https://upload.wikimedia.org/wikipedia/commons/3/39/GodfreyKneller-IsaacNewton-1689.jpg)
background-position: 95% 10%  
background-size: 20%

# Newton’s Method

Our curent "guess" for a: `\(x_0\)` ➡️ `\(f'(x_0) = \frac{f(x_0) - y}{x_0 - x}\)`   
➡️ `\(f'(x_0) = \frac{f(x_0) - 0}{x_0 - x_1}\)` ➡️ `\(x_1= x_0 - \frac{f(x_0)}{f'(x_0)}\)` ➡️ `\(x_2= x_1 - \frac{f(x_1)}{f'(x_1)}\)`

.center[
<img src="newton_method.png" width="60%" height="60%" />
] 
---
# Newton’s Method

.center[
![](https://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif)]

---
# Newton’s Method

### In general: 

## `\(x_{n+1}= x_n - \frac{f(x_n)}{f'(x_n)}\)`

---
# Newton’s Method

###  Computational example:

### Find roots for `\(x^3-1\)`

```r
F1 <- function(x){
 return(c(x^3-1, 3*(x^2))) 
} 
```
---
# Newton’s Method

###  Computational example:

```r
newtonraphson <- 
 function(ftn, x0, tol = 1e-9, max.iter) {
 # initialize 
 x <- x0
 fx <- ftn(x)
 iter <- 0
 
 # continue iterating until stopping conditions are met 
 while((abs(fx[1]) > tol) && (iter < max.iter)) {
 x <- x - fx[1]/fx[2]
 fx <- ftn(x)
 iter <- iter + 1
 cat("At iteration", iter, "value of x is:", x, "\n")
 }
 
 # output depends on the success of the algorithm 
 if (abs(fx[1]) > tol){
 cat("Algorithm failed to converge\n")
 return(data.frame(x0, root = NA, iter = NA))
 } else {
 cat("Algorithm converged\n")
 return(data.frame(x0, root = x, iter))
 }
```

---
background-image: url(https://upload.wikimedia.org/wikipedia/commons/thumb/0/06/CIRCLE_LINES.svg/594px-CIRCLE_LINES.svg.png)
background-position: 95% 8%  
background-size: 20%   
# Secant Method

- Do not need to compute a derivative   
- Need to provide two initial guesses

---
background-image: url(https://upload.wikimedia.org/wikipedia/commons/thumb/0/06/CIRCLE_LINES.svg/594px-CIRCLE_LINES.svg.png)
background-position: 95% 8%  
background-size: 20%   
# Secant Method        
Two initial "guesses" `\(x_0\)` and `\(x_1\)`, assuming `\(x_0\)` is the older one   
➡️ `\(\frac{y-f(x_1)}{x-x_1} = \frac{f(x_0) - f(x_1)}{x_0 - x_1}\)`   
➡️ so `\(x_2\)` can be found from `\(\frac{0-f(x_1)}{x_2-x_1} = \frac{f(x_0) - f(x_1)}{x_0 - x_1}\)` ➡️  `\(x_2 = x_1 - f(x_1)\frac{x_0 - x_1}{f(x_0) - f(x_1)}\)`

---

# Secant Method 

### In general: 

## `\(x_{n+1} = x_n - f(x_n)\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}\)`

---

# Secant Method

###  Computational example:

```r
secant <- function(ftn, x0, x1, tol = 1e-9, max.iter) {
 # initialize 
 x_n0 <- x0
 x_n1 <- x1
 ftn_n0 <- ftn(x_n0)
 ftn_n1 <- ftn(x_n1)
 iter <- 0
 
 # continue iterating until stopping conditions are met 
 while((abs(ftn_n1) > tol) && (iter < max.iter)) {
 x_n2 <- x_n1 - ftn_n1*(x_n1 - x_n0)/(ftn_n1 - ftn_n0)
 x_n0 <- x_n1
 ftn_n0 <- ftn(x_n0)
 x_n1 <- x_n2
 ftn_n1 <- ftn(x_n1)
 iter <- iter + 1
 cat("At iteration", iter, "value of x is:", x_n1, "\n")
 }
 
 return(c(x_n1, iter))
}
```

---

# Fractals

### Root finding functions can also be applied to find roots for **complex functions**, which are functions of **complex numbers**.

<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/a/af/Complex_number_illustration.svg/1024px-Complex_number_illustration.svg.png" width="50%" height="50%" />
---
# Fractals

### Each root has a **basin of attraction** in the complex plane, which is a set of all **initial guesses** that cause the method to converge to that particular root.

### These sets can be mapped into images. The boundaries of the basins of attraction are called **fractals**.

---
# Fractals

## Julia set

---
# Fractals

## Mandelbrot set

<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/2/21/Mandel_zoom_00_mandelbrot_set.jpg/1920px-Mandel_zoom_00_mandelbrot_set.jpg" width="50%" height="50%" />
---
# Fractals

## Burning ship fractals

<img src="https://upload.wikimedia.org/wikipedia/commons/9/9a/Burning_Ship_Fractal_Zoom-out_64.gif" width="50%" height="50%" />
---
# Fractals

## Newton's fractals

<img src="https://upload.wikimedia.org/wikipedia/commons/9/9b/FRACT008.png" width="50%" height="50%" />
---
# Summary

### 1. Root-finding algorithm: `\(h(x) = f(x) – g(x) = 0\)` 📐

### 2. Most common: Newton’s Method and Secant method   🍎

### 3. Complex numbers and fractals 👩‍🎨

---
# Get Started - Creating your fractals 
 
.center[
![](https://media1.tenor.com/images/58863bf1b0b453857aa0a3e390eac3c6/tenor.gif?itemid=5675661) 
]

---
# Recommended readings  
-  Introduction to scientific programming and simulation using r, by Andrew P. Robinson, Owen Jones, and Robert Maillardet [[link](http://www.tf.uns.ac.rs/~omorr/radovan_omorjan_003_prII/r-examples/spuRs/spuRs-2ed.pdf)]

- Newton fractal wiki page  [[link](https://en.wikipedia.org/wiki/Newton_fractal)]

---

class: center, middle,inverse  
# Thanks for attending!

😄

### Keep in touch [twitter @kerenxuepi](https://twitter.com/kerenxuepi) 
 
Slides created via **RLadies** **xaringan** slide theme