daviddoji/CalculusWithJuliaNotes.jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'main' of https://github.com/jverzani/CalculusWithJuliaNotes.jl

Browse Source
This commit is contained in:
jverzani 2023-05-05 11:28:00 -04:00
commit aeb3364add
1 changed files with 12 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
quarto/derivatives/more_zeros.qmd
View File

@ -54,7 +54,7 @@ One way to view Newton's method is through the inverse of $f$ (assuming it exist
If $f$ has a simple zero at $\alpha$ and is locally invertible (that is some $f^{-1}$ exists) then the update step for Newton's method can be identified with:
* fitting a polynomial to the local inverse function of $f$ going through through the point $(f(x_0),x_0)$,
* fitting a polynomial to the local inverse function of $f$ going through the point $(f(x_0),x_0)$,
* and matching the slope of $f$ at the same point.
@ -67,7 +67,7 @@ A similar viewpoint can be used to create derivative-free methods.
For example, the [secant method](https://en.wikipedia.org/wiki/Secant_method) can be seen as the result of fitting a degree-$1$ polynomial approximation for $f^{-1}$ through two points $(f(x_0),x_0)$ and $(f(x_1), x_1)$.
Again, expressing this approximation as $g(y) = h_0 + h_1(y-f(x_1))$ leads to $g(f(x_1)) = x_1 = h_0$. Substituting $f(x_0)$ gives $g(f(x_0)) = x_0 = x_1 + h_1(f(x_0)-f(x_1))$. Solving for $h_1$ leads to $h_1=(x_1-x_0)/(f(x_1)-f(x_0))$. Then $x_2 = g(0) = x_1 + (x_1-x_0)/(f(x_1)-f(x_0)) \cdot f(x_1)$. This is the first step of the secant method:
Again, expressing this approximation as $g(y) = h_0 + h_1(y-f(x_1))$ leads to $g(f(x_1)) = x_1 = h_0$. Substituting $f(x_0)$ gives $g(f(x_0)) = x_0 = x_1 + h_1(f(x_0)-f(x_1))$. Solving for $h_1$ leads to $h_1=(x_1-x_0)/(f(x_1)-f(x_0))$. Then $x_2 = g(0) = x_1 - (x_1-x_0)/(f(x_1)-f(x_0)) \cdot f(x_1)$. This is the first step of the secant method:
$$
@ -129,7 +129,7 @@ Let $\epsilon_{n+1} = x_{n+1}-\alpha$, where $\alpha$ is assumed to be the *simp
\begin{align*}
\epsilon_{n+1} &\approx \frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})} \frac{(1/2)f''(\alpha)(e_n-e_{n-1})}{x_n-x_{n-1}} \epsilon_n \epsilon_{n-1}\\
\epsilon_{n+1} &\approx \frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})} \frac{(1/2)f''(\alpha)(\epsilon_n-\epsilon_{n-1})}{x_n-x_{n-1}} \epsilon_n \epsilon_{n-1}\\
& \approx \frac{f''(\alpha)}{2f'(\alpha)} \epsilon_n \epsilon_{n-1}\\
&= C \epsilon_n \epsilon_{n-1}.
\end{align*}
@ -142,7 +142,7 @@ Assuming the error term has the form $\epsilon_{n+1} = A|\epsilon_n|^\phi$ and s
$$
\frac{A^{1-1/\phi}}{C} = |\epsilon_n|^{1 - \phi +1/\phi}.
\frac{A^{1+1/\phi}}{C} = |\epsilon_n|^{1 - \phi +1/\phi}.
$$
The left side being a constant suggests $\phi$ solves: $1 - \phi + 1/\phi = 0$ or $\phi^2 -\phi - 1 = 0$. The solution is the golden ratio, $(1 + \sqrt{5})/2 \approx 1.618\dots$.
@ -207,7 +207,7 @@ Though the above can be simplified quite a bit when computed by hand, here we si
(`SymPy`'s `lambdify` function, by default, picks the order of its argument lexicographically, in this case they will be the `f` values then the `x` values.)
An inverse quadratic step is utilized by Brent's method, as possible, to yield a rapidly convergent bracketing algorithm implemented as a default zero finder in many software languages. `Julia`'s `Roots` package implements the method in `Roots.Brent()`. An inverse cubic interpolation is utilized by [Alefeld, Potra, and Shi](https://dl.acm.org/doi/10.1145/210089.210111) which gives an asymptotically even more rapidly convergent algorithm then Brent's (implemented in `Roots.AlefeldPotraShi()` and also `Roots.A42()`). This is used as a finishing step in many cases by the default hybrid `Order0()` method of `find_zero`.
An inverse quadratic step is utilized by Brent's method, as possible, to yield a rapidly convergent bracketing algorithm implemented as a default zero finder in many software languages. `Julia`'s `Roots` package implements the method in `Roots.Brent()`. An inverse cubic interpolation is utilized by [Alefeld, Potra, and Shi](https://dl.acm.org/doi/10.1145/210089.210111) which gives an asymptotically even more rapidly convergent algorithm than Brent's (implemented in `Roots.AlefeldPotraShi()` and also `Roots.A42()`). This is used as a finishing step in many cases by the default hybrid `Order0()` method of `find_zero`.
In a bracketing algorithm, the next step should reduce the size of the bracket, so the next iterate should be inside the current bracket. However, quadratic convergence does not guarantee this to happen. As such, sometimes a subsitute method must be chosen.
@ -371,10 +371,10 @@ However, there may be no floating point value where $f$ is exactly `0.0` so chec
First if `f(x_n)` is `0.0` then it makes sense to call `x_n` an *exact zero* of $f$, even though this may hold even if `x_n`, a floating point value, is not mathematically an *exact* zero of $f$. (Consider `f(x) = x^2 - 2x + 1`. Mathematically, this is identical to `g(x) = (x-1)^2`, but `f(1 + eps())` is zero, while `g(1+eps())` is `4.930380657631324e-32`.
However, there may never be a value with `f(x_n)` exactly `0.0`. (The value of `sin(pi)` is not zero, for example, as `pi` is an approximation to $\pi$, as well the `sin` of values adjacent to `float(pi)` do not produce `0.0` exactly.)
However, there may never be a value with `f(x_n)` exactly `0.0`. (The value of `sin(1pi)` is not zero, for example, as `1pi` is an approximation to $\pi$, as well the `sin` of values adjacent to `float(pi)` do not produce `0.0` exactly.)
Suppose `x_n` is the closest floating number to $\alpha$, the zero. Then the relative rounding error, $($ `x_n` $- \alpha)/\alpha$, will be a value $(1 + \delta)$ with $\delta$ less than `eps()`.
Suppose `x_n` is the closest floating number to $\alpha$, the zero. Then the relative rounding error, $($ `x_n` $- \alpha)/\alpha$, will be a value $\delta$ with $\delta$ less than `eps()`.
How far then can `f(x_n)` be from $0 = f(\alpha)$?
@ -491,7 +491,7 @@ numericq(10, 1)
###### Question
Let `f(x) = exp(x) - x^4` and a starting bracket be `x0 = [8.9]`. Then calling `find_zero(f,x0, verbose=true)` will show that 49 steps are needed for exact bisection to converge. What about with the `Roots.Brent()` algorithm, which uses inverse quadratic steps when it can?
Let `f(x) = exp(x) - x^4` and a starting bracket be `x0 = [8, 9]`. Then calling `find_zero(f,x0, verbose=true)` will show that 48 steps are needed for exact bisection to converge. What about with the `Roots.Brent()` algorithm, which uses inverse quadratic steps when it can?
It takes how many steps?
@ -509,10 +509,10 @@ The `Roots.A42()` method uses inverse cubic interpolation, as possible, how many
```{julia}
#| hold: true
#| echo: false
numericq(3, 1)
numericq(7, 1)
```
The large difference is due to how the tolerances are set within `Roots`. The `Brent method gets pretty close in a few steps, but takes a much longer time to get close enough for the default tolerances
The large difference is due to how the tolerances are set within `Roots`. The Brent method gets pretty close in a few steps, but takes a much longer time to get close enough for the default tolerances.
###### Question
@ -526,7 +526,7 @@ Consider this crazy function defined by:
f(x) = cos(100*x)-4*erf(30*x-10)
```
(The `erf` function is the (error function](https://en.wikipedia.org/wiki/Error_function) and is in the `SpecialFunctions` package loaded with `CalculusWithJulia`.)
(The `erf` function is the [error function](https://en.wikipedia.org/wiki/Error_function) and is in the `SpecialFunctions` package loaded with `CalculusWithJulia`.)
Make a plot over the interval $[-3,3]$ to see why it is called "crazy".
@ -568,7 +568,7 @@ answ = 4
radioq(choices, answ, keep_order=true)
```
Does `find_zero` find a zero to this function starting from $1$?
Does `find_zero` find a zero to this function starting from $0.175$?
```{julia}