diff --git a/quarto/_quarto.yml b/quarto/_quarto.yml index 40dada1..d44c3cd 100644 --- a/quarto/_quarto.yml +++ b/quarto/_quarto.yml @@ -138,7 +138,9 @@ website: format: html: - theme: lux # spacelab # lux # sketchy # cosmo # https://quarto.org/docs/output-formats/html-themes.html + theme: + light: lux # spacelab # lux # sketchy # cosmo # https://quarto.org/docs/output-formats/html-themes.html + dark: darkly number-depth: 3 toc-depth: 3 link-external-newwindow: true diff --git a/quarto/derivatives/curve_sketching.qmd b/quarto/derivatives/curve_sketching.qmd index fd79296..48f8b1f 100644 --- a/quarto/derivatives/curve_sketching.qmd +++ b/quarto/derivatives/curve_sketching.qmd @@ -12,7 +12,7 @@ using Plots plotly() using SymPy using Roots -using Polynomials # some name clash with SymPy +import Polynomials: variable, Polynomial, RationalFunction # avoid name clash ``` @@ -162,8 +162,8 @@ This sort of analysis can be automated. The plot "recipe" for polynomials from t ```{julia} -xβ‚š = variable(Polynomial) -plot(f(xβ‚š)) # f(xβ‚š) of Polynomial type +𝐱 = variable(Polynomial) +plot(f(𝐱)) # f(𝐱) of Polynomial type ``` ##### Example @@ -219,8 +219,8 @@ Again, this sort of analysis can be systematized. The rational function type in ```{julia} -xα΅£ = variable(RationalFunction) -plot(f(xα΅£)) # f(x) of RationalFunction type +𝐱 = variable(RationalFunction) +plot(f(𝐱)) # f(x) of RationalFunction type ``` ##### Example diff --git a/quarto/derivatives/derivatives.qmd b/quarto/derivatives/derivatives.qmd index 1160ffd..d60fbc1 100644 --- a/quarto/derivatives/derivatives.qmd +++ b/quarto/derivatives/derivatives.qmd @@ -730,13 +730,13 @@ $$ \frac{f(g(x+h)) - f(g(x))}{h} = \frac{f(g(x+h)) - f(g(x))}{g(x+h) - g(x)} \cdot \frac{g(x+h) - g(x)}{h}. $$ -The left hand side will converge to the derivative of $u(x)$ or $[f(g(x))]'$. +The left-hand side will converge to the derivative of $u(x)$ or $[f(g(x))]'$. -The right most part of the right side would have a limit $g'(x)$, were we to let $h$ go to $0$. +The right-most part of the right-hand side would have a limit $g'(x)$, were we to let $h$ go to $0$. -It isn't obvious, but the left part of the right side has the limit $f'(g(x))$. This would be clear if *only* $g(x+h) = g(x) + h$, for then the expression would be exactly the limit expression with $c=g(x)$. But, alas, except to some hopeful students and some special cases, it is definitely not the case in general that $g(x+h) = g(x) + h$ - that right parentheses actually means something. However, it is *nearly* the case that $g(x+h) = g(x) + kh$ for some $k$ and this can be used to formulate a proof (one of the two detailed [here](http://en.wikipedia.org/wiki/Chain_rule#Proofs) and [here](http://kruel.co/math/chainrule.pdf)). +It isn't obvious, but the left part of the right-hand side has the limit $f'(g(x))$. This would be clear if *only* $g(x+h) = g(x) + h$, for then the expression would be exactly the limit expression with $c=g(x)$. But, alas, except to some hopeful students and some special cases, it is definitely not the case in general that $g(x+h) = g(x) + h$ - that right parentheses actually means something. However, it is *nearly* the case that $g(x+h) = g(x) + kh$ for some $k$ and this can be used to formulate a proof (one of the two detailed [here](http://en.wikipedia.org/wiki/Chain_rule#Proofs) and [here](http://kruel.co/math/chainrule.pdf)). Combined, we would end up with: diff --git a/quarto/derivatives/first_second_derivatives.qmd b/quarto/derivatives/first_second_derivatives.qmd index 9f8c749..9baad5a 100644 --- a/quarto/derivatives/first_second_derivatives.qmd +++ b/quarto/derivatives/first_second_derivatives.qmd @@ -162,9 +162,10 @@ This question can be answered by considering the first derivative. > *The first derivative test*: If $c$ is a critical point for $f(x)$ and *if* $f'(x)$ changes sign at $x=c$, then $f(c)$ will be either a relative maximum or a relative minimum. > -> * It will be a relative maximum if the derivative changes sign from $+$ to $-$. -> * It will be a relative minimum if the derivative changes sign from $-$ to $+$. -> * If $f'(x)$ does not change sign at $c$, then $f(c)$ is *not* a relative maximum or minimum. +> * $f$ will have a relative maximum at $c$ if the derivative changes sign from $+$ to $-$. +> * $f$ will have a relative minimum at $c$ if the derivative changes sign from $-$ to $+$. +> +> Further, If $f'(x)$ does *not* change sign at $c$, then $f$ will *not* have a relative maximum or minimum at $c$. @@ -185,7 +186,7 @@ f(x) = exp(-abs(x)) * cos(pi * x) plotif(f, f', -3, 3) ``` -We can see the first derivative test in action: at the peaks and valleys – the relative extrema – the color changes. This is because $f'$ is changing sign as the function changes from increasing to decreasing or vice versa. +We can see the first derivative test in action: at the peaks and valleys – the relative extrema – the highlighting changes. This is because $f'$ is changing sign as the function changes from increasing to decreasing or vice versa. This function has a critical point at $0$, as can be seen. It corresponds to a point where the derivative does not exist. It is still identified through `find_zeros`, which picks up zeros and in case of discontinuous functions, like `f'`, zero crossings: @@ -493,8 +494,8 @@ Concave up functions are "opening" up, and often clearly $U$-shaped, though that > The **second derivative test**: If $c$ is a critical point of $f(x)$ with $f''(c)$ existing in a neighborhood of $c$, then > -> * The value $f(c)$ will be a relative maximum if $f''(c) > 0$, -> * The value $f(c)$ will be a relative minimum if $f''(c) < 0$, and +> * $f$ will have a relative maximum at the critical point $c$ if $f''(c) > 0$, +> * $f$ will have a relative minimum at the critical point $c$ if $f''(c) < 0$, and > * *if* $f''(c) = 0$ the test is *inconclusive*. diff --git a/quarto/derivatives/linearization.qmd b/quarto/derivatives/linearization.qmd index 1ab5e63..ef877cb 100644 --- a/quarto/derivatives/linearization.qmd +++ b/quarto/derivatives/linearization.qmd @@ -40,7 +40,11 @@ tangent(f, c) = x -> f(c) + f'(c) * (x - c) (Recall, the `->` indicates that an anonymous function is being generated.) -This function along with the `f'` notation for automatic derivatives is defined in the `CalculusWithJulia` package. +This function along with the `f'` notation for automatic derivatives is defined in the `CalculusWithJulia` package, though a bit differently: + +```{julia} +tangent(sin, pi/4) +``` We make some graphs with tangent lines: diff --git a/quarto/derivatives/mean_value_theorem.qmd b/quarto/derivatives/mean_value_theorem.qmd index 318f855..64222c9 100644 --- a/quarto/derivatives/mean_value_theorem.qmd +++ b/quarto/derivatives/mean_value_theorem.qmd @@ -306,7 +306,7 @@ annotate!([(0,jβ‚€,text("a", :bottom)), ## The mean value theorem -We are driving south and in one hour cover 70 miles. If the speed limit is 65 miles per hour, were we ever speeding? We'll we averaged more than the speed limit so we know the answer is yes, but why? Speeding would mean our instantaneous speed was more than the speed limit, yet we only know for sure our *average* speed was more than the speed limit. The mean value tells us that if some conditions are met, then at some point (possibly more than one) we must have that our instantaneous speed is equal to our average speed. +We are driving south and in one hour cover 70 miles. If the speed limit is 65 miles per hour, were we ever speeding? Well we averaged more than the speed limit so we know the answer is yes, but why? Speeding would mean our instantaneous speed was more than the speed limit, yet we only know for sure our *average* speed was more than the speed limit. The mean value tells us that if some conditions are met, then at some point (possibly more than one) we must have that our instantaneous speed is equal to our average speed. The mean value theorem is a direct generalization of Rolle's theorem. @@ -404,7 +404,7 @@ board.create('tangent', [r], {strokeColor:'#ff0000'}); line = board.create('line',[p[0],p[1]],{strokeColor:'#ff0000',dash:1}); ``` -This interactive example can also be found at [jsxgraph](http://jsxgraph.uni-bayreuth.de/wiki/index.php?title=Mean_Value_Theorem). It shows a cubic polynomial fit to the $4$ adjustable points labeled A through D. The secant line is drawn between points A and B with a dashed line. A tangent line – with the same slope as the secant line – is identified at a point $(\alpha, f(\alpha))$ where $\alpha$ is between the points A and B. That this can always be done is a conseuqence of the mean value theorem. +This interactive example can also be found at [jsxgraph](http://jsxgraph.uni-bayreuth.de/wiki/index.php?title=Mean_Value_Theorem). It shows a cubic polynomial fit to the $4$ adjustable points labeled A through D. The secant line is drawn between points A and B with a dashed line. A tangent line – with the same slope as the secant line – is identified at a point $(\alpha, f(\alpha))$ where $\alpha$ is between the points A and B. That this can always be done is a consequence of the mean value theorem. ##### Example diff --git a/quarto/derivatives/newtons_method.qmd b/quarto/derivatives/newtons_method.qmd index 98631cd..eb0181f 100644 --- a/quarto/derivatives/newtons_method.qmd +++ b/quarto/derivatives/newtons_method.qmd @@ -864,7 +864,7 @@ The calculation that produces the quadratic convergence now becomes: $$ \begin{align*} -x_{i+1} - \alpha &= (x_i - \alpha) - \frac{1}{k}(x_i-\alpha - \frac{f''(\xi)}{2f'(x_i)}(x_i-\alpha)^2) +x_{i+1} - \alpha &= (x_i - \alpha) - \frac{1}{k}(x_i-\alpha - \frac{f''(\xi)}{2f'(x_i)}(x_i-\alpha)^2) \\ &= \frac{k-1}{k} (x_i-\alpha) + \frac{f''(\xi)}{2kf'(x_i)}(x_i-\alpha)^2. \end{align*} $$ diff --git a/quarto/derivatives/numeric_derivatives.qmd b/quarto/derivatives/numeric_derivatives.qmd index 8606978..b144cb7 100644 --- a/quarto/derivatives/numeric_derivatives.qmd +++ b/quarto/derivatives/numeric_derivatives.qmd @@ -195,11 +195,10 @@ We pass in the function object, `k''`, and not the evaluated function. ## Recap on derivatives in Julia -A quick summary for finding derivatives in `Julia`, as there are $3$ different manners: - +A quick summary of the $3$ different ways for finding derivatives in `Julia` presented in these notes: * Symbolic derivatives are found using `diff` from `SymPy` - * Automatic derivatives are found using the notation `f'` using `ForwardDiff.derivative` + * Automatic derivatives are found using the notation `f'` which utilizes `ForwardDiff.derivative` * approximate derivatives at a point, `c`, for a given `h` are found with `(f(c+h)-f(c))/h`. diff --git a/quarto/derivatives/optimization.qmd b/quarto/derivatives/optimization.qmd index 11aed8f..cf13234 100644 --- a/quarto/derivatives/optimization.qmd +++ b/quarto/derivatives/optimization.qmd @@ -38,7 +38,7 @@ The main tool is the extreme value theorem of Bolzano and Fermat's theorem about Though not all of our problems lend themselves to a description of a continuous function on a closed interval, if they do, we have an algorithmic prescription to find the absolute extrema of a function: -1. Find the critical points. For this we can use a root-finding algorithm like `find_zero`. +1. Find the critical points. For this we can use a root-finding algorithm like `find_zero` or `find_zeros`. 2. Evaluate the function values at the critical points and at the end points. 3. Identify the largest and smallest values. @@ -65,7 +65,9 @@ gr() function perimeter_area_graphic_graph(n) h = 1 + 2n w = 10-h - plt = plot([0,0,w,w,0], [0,h,h,0,0], legend=false, size=fig_size, + plt = plot(Shape([0,0,w,w,0], [0,h,h,0,0]), legend=false, + fill=(:gray, 0.1), + size=fig_size, xlim=(0,10), ylim=(0,10)) scatter!(plt, [w], [h], color=:orange, markersize=5) annotate!(plt, [(w/2, h/2, "Area=$(round(w*h,digits=1))")]) @@ -153,12 +155,16 @@ Now we can solve graphically as before, or numerically, such as here where we se find_zeros(A', 0, 10) # find_zeros in `Roots`, ``` -(As a reminder, the notation `A'` is defined in `CalculusWithJulia` using the `derivative` function from the `ForwardDiff` package.) - :::{.callout-note} ## Note -Look at the last definition of `A`. The function `A` appears on both sides, though on the left side with one argument and on the right with two. These are two "methods" of a *generic* function, `A`. `Julia` allows multiple definitions for the same name as long as the arguments (their number and type) can disambiguate which to use. In this instance, when one argument is passed in then the last definition is used (`A(b,h(b))`), whereas if two are passed in, then the method that multiplies both arguments is used. The advantage of multiple dispatch is illustrated: the same concept - area - has one function name, though there may be different ways to compute the area, so there is more than one implementation. +Look at the last definition of `A`. The function `A` appears on both sides, though on the left side with one argument and on the right with two. These are two "methods" of a *generic* function, `A`. `Julia` allows multiple definitions (methods) for the same function name (a generic function) as long as the arguments (their number and type) can disambiguate which to use. In this instance, when one argument is passed in then the last definition is used (`A(b, h(b))`), whereas if two are passed in, then the method that multiplies both arguments is used. The advantage of multiple dispatch is illustrated: the same concept - area - has one function name, though there may be different ways to compute the area, so there is more than one implementation. + +::: + +:::{.callout-note} +## Alternatively ... +If using composition for substitution is too complicated, the step of creating `A(b,h)` can be skipped by just using `h(b)` when defining `A`: `A(b) = b * h(b)`. ::: @@ -203,7 +209,7 @@ In `Julia` this is ```{julia} -Aα΅£(x, y) = x*y + pi*(x/2)^2 / 2 +A(x, y) = x*y + pi*(x/2)^2 / 2 ``` The perimeter consists of $3$ sides of the rectangle and the perimeter of half a circle ($\pi r$, with $r=x/2$): @@ -224,7 +230,7 @@ And then we substitute in `y(x)` for `y` in the area formula through: ```{julia} -Aα΅£(x) = Aα΅£(x, y(x)) +A(x) = A(x, y(x)) ``` Of course both $x$ and $y$ are non-negative. The latter forces $x$ to be no more than $x=20/(1+\pi/2)$. @@ -237,7 +243,7 @@ We begin by simply graphing and estimating the values of the maximum and where i ```{julia} -plot(Aα΅£, 0, 20/(1+pi/2)) +plot(A, 0, 20/(1+pi/2)) ``` The naked eye sees that maximum value is somewhere around $27$ and occurs at $x\approx 5.6$. Clearly from the graph, we know the maximum value happens at the critical point and there is only one such critical point. @@ -247,7 +253,7 @@ As reading the maximum from the graph is more difficult than reading a $0$ of a ```{julia} -plot(Aα΅£', 5.5, 5.7) +plot(A', 5.5, 5.7) ``` We confirm that the critical point is around $5.6$. @@ -259,18 +265,21 @@ We confirm that the critical point is around $5.6$. Rather than zoom in graphically, we now use a root-finding algorithm, to find a more precise value of the zero of $A'$. We know that the maximum will occur at a critical point, a zero of the derivative. The `find_zero` function from the `Roots` package provides a non-linear root-finding algorithm based on the bisection method. The only thing to keep track of is that solving $f'(x) = 0$ means we use the derivative and not the original function. -We see from the graph that $[0, 20/(1+\pi/2)]$ will provide a bracket, as there is only one relative maximum: +In a few sections we will see how to use `find_zero` to find a zero *near* an initial guess, but for now we will use `find_zero` to find a zero *between* two values which form a bracketing interval. + + +From the graph, it is clear that $[0, 20/(1+\pi/2)]$ will provide a bracket for the derivative, as there is only one relative maximum: ```{julia} -xβ€² = find_zero(Aα΅£', (0, 20/(1+pi/2))) +xβ€² = find_zero(A', (0, 20/(1+pi/2))) ``` This value is the lone critical point, and in this case gives the position of the value that will maximize the function. The value and maximum area are then given by: ```{julia} -(xβ€², Aα΅£(xβ€²)) +(xβ€², A(xβ€²)) ``` (Compare this answer to the previous, is the square the figure of greatest area for a fixed perimeter, or just the figure amongst all rectangles? See [Isoperimetric inequality](https://en.wikipedia.org/wiki/Isoperimetric_inequality) for an answer.) @@ -279,11 +288,11 @@ This value is the lone critical point, and in this case gives the position of th ### Using `argmax` to identify where a function is maximized -This value that maximizes a function is sometimes referred to as the *argmax*, or argument which maximizes the function. In `Julia` the `argmax(f,domain)` function is defined to "Return a value $x$ in the domain of $f$ for which $f(x)$ is maximized. If there are multiple maximal values for $f(x)$ then the first one will be found." The domain is some iterable collection. In the mathematical world this would be an interval $[a,b]$, but on the computer it is an approximation, such as is returned by `range` below. Without having to take a derivative, as above, but sacrificing some accuracy, the task of identifying `x` for where `A` is maximum, could be done with +This value that maximizes a function is sometimes referred to as the *argmax*, or argument which maximizes the function. In `Julia` the `argmax(f, domain)` function is defined to "Return a value $x$ in the domain of $f$ for which $f(x)$ is maximized. If there are multiple maximal values for $f(x)$ then the first one will be found." The domain is some iterable collection. In the mathematical world this would be an interval $[a,b]$, but on the computer it is an approximation, such as is returned by `range` below. Without having to take a derivative, as above, but sacrificing some accuracy, the task of identifying `x` for where `A` is maximum, could be done with ```{julia} -argmax(Aα΅£, range(0, 20/(1+pi/2), length=10000)) +argmax(A, range(0, 20/(1+pi/2), length=10000)) ``` #### A symbolic approach @@ -293,20 +302,20 @@ We could also do the above problem symbolically with the aid of `SymPy`. Here ar ```{julia} -@syms π’˜::real 𝒉::real +@syms 𝐰::real 𝐑::real -𝑨₀ = π’˜ * 𝒉 + pi * (π’˜/2)^2 / 2 -𝑷erim = 2*𝒉 + π’˜ + pi * π’˜/2 -𝒉₀ = solve(𝑷erim - 20, 𝒉)[1] -𝑨₁ = 𝑨₀(𝒉 => 𝒉₀) -π’˜β‚€ = solve(diff(𝑨₁,π’˜), π’˜)[1] +𝐀₀ = 𝐰 * 𝐑 + pi * (𝐰/2)^2 / 2 +𝐏erim = 2*𝐑 + 𝐰 + pi * 𝐰/2 +𝐑₀ = solve(𝐏erim - 20, 𝐑)[1] +𝐀₁ = 𝐀₀(𝐑 => 𝐑₀) +𝐰₀ = solve(diff(𝐀₁,𝐰), 𝐰)[1] ``` -We know that `π’˜β‚€` is the maximum in this example from our previous work. We shall see soon, that just knowing that the second derivative is negative at `π’˜β‚€` would suffice to know this. Here we check that condition: +We know that `𝐰₀` is the maximum in this example from our previous work. We shall see soon, that just knowing that the second derivative is negative at `𝐰₀` would suffice to know this. Here we check that condition: ```{julia} -diff(𝑨₁, π’˜, π’˜)(π’˜ => π’˜β‚€) +diff(𝐀₁, 𝐰, 𝐰)(𝐰 => 𝐰₀) ``` As an aside, compare the steps involved above for a symbolic solution to those of previous work for a numeric solution: @@ -314,9 +323,9 @@ As an aside, compare the steps involved above for a symbolic solution to those o ```{julia} #| hold: true -Aα΅£(w, h) = w*h + pi*(w/2)^2 / 2 -h(w) = (20 - w - pi * w/2) / 2 -Aα΅£(w) = A(w, h(w)) +A(w, h) = w * h + pi*(w/2)^2 / 2 +h(w) = (20 - w - pi * w/2) / 2 +A(w) = A(w, h(w)) find_zero(A', (0, 20/(1+pi/2))) # 40 / (pi + 4) ``` @@ -496,7 +505,7 @@ It appears to be around $8.3$. We now use `find_zero` to refine our guess at the Ξ± = find_zero(T', (7, 9)) ``` -Okay, got it. Around$8.23$. So is our minimum time +Okay, got it. Around $8.23$. So is our minimum time. ```{julia} @@ -527,9 +536,9 @@ caption = L""" Image number $43$ from l'Hospital's calculus book (the first). A traveler leaving location $C$ to go to location $F$ must cross two regions separated by the straight line $AEB$. We suppose that in the -region on the side of $C$, he covers distance $a$ in time $c$, and +region on the side of $C$, they cover distance $a$ in time $c$, and that on the other, on the side of $F$, distance $b$ in the same time -$c$. We ask through which point $E$ on the line $AEB$ he should pass, +$c$. We ask through which point $E$ on the line $AEB$ they should pass, so as to take the least possible time to get from $C$ to $F$? (From http://www.ams.org/samplings/feature-column/fc-2016-05.) @@ -675,7 +684,9 @@ Here traveling directly to the point $(L,0)$ is fastest. Though travel is slower Maximize the function $xe^{-(1/2) x^2}$ over the interval $[0, \infty)$. -Here the extreme value theorem doesn't technically apply, as we don't have a closed interval. However, **if** we can eliminate the endpoints as candidates, then we should be able to convince ourselves the maximum must occur at a critical point of $f(x)$. (If not, then convince yourself for all sufficiently large $M$ the maximum over $[0,M]$ occurs at a critical point, not an endpoint. Then let $M$ go to infinity. In general, for an optimization problem of a continuous function on the interval $(a,b)$ if the right limit at $a$ and left limit at $b$ can be ruled out as candidates, the optimal value must occur at a critical point.) +Here the extreme value theorem doesn't technically apply, as we don't have a closed interval. However, **if** we can eliminate the endpoints as candidates, then we should be able to convince ourselves the maximum must occur at a critical point of $f(x)$. (If not, then convince yourself for all sufficiently large $M$ the maximum over $[0,M]$ occurs at a critical point, not an endpoint. Then let $M$ go to infinity. + +In general, for an optimization problem of a continuous function on the interval $(a,b)$ if the right limit at $a$ and left limit at $b$ can be ruled out as candidates, the optimal value must occur at a critical point.) So to approach this problem we first graph it over a wide interval. @@ -762,14 +773,14 @@ The minimum looks to be around $4$cm and is clearly between $2$cm and $6$cm. We ```{julia} -rₛₐ = find_zero(SA', (2, 6)) +r = find_zero(SA', (2, 6)) ``` Okay, $3.837...$ is our answer for $r$. Use this to get $h$: ```{julia} -canheight(rₛₐ) +canheight(r) ``` This produces a can which is about square in profile. This is not how most cans look though. Perhaps our model is too simple, or the cans are optimized for some other purpose than minimizing materials. @@ -1417,7 +1428,7 @@ radioq(choices, 1) ###### Question -In [Hall](https://www.maa.org/sites/default/files/hall03010308158.pdf) we can find a dozen optimization problems related to the following figure of the parabola $y=x^2$ a point $P=(a,a^2)$, $a > 0$, and its normal line. We will do two. +In [Hall](https://www.maa.org/sites/default/files/hall03010308158.pdf) we can find a dozen optimization problems related to the following figure of the parabola $y=x^2$ a point $P=(a,a^2)$, $a > 0$, and its normal line. We will do two here, the full set will be looked at later. ```{julia} diff --git a/quarto/derivatives/symbolic_derivatives.qmd b/quarto/derivatives/symbolic_derivatives.qmd index a20fe19..1308da3 100644 --- a/quarto/derivatives/symbolic_derivatives.qmd +++ b/quarto/derivatives/symbolic_derivatives.qmd @@ -3,7 +3,7 @@ {{< include ../_common_code.qmd >}} -This section uses this add-on package: +This section uses the `TermInterface` add-on package. ```{julia} @@ -54,16 +54,7 @@ arguments(:(-x)), arguments(:(pi^2)), arguments(:(1 + x + x^2)) (The last one may be surprising, but all three arguments are passed to the `+` function.) -Here we define a function to decide the *arity* of an expression based on the number of arguments it is called with: - - -```{julia} -function arity(ex) - n = length(arguments(ex)) - n == 1 ? Val(:unary) : - n == 2 ? Val(:binary) : Val(:nary) -end -``` +`TermInterface` has an `arity` function defined by `length(arguments(ex))` that will be used for dispatch below. Differentiation must distinguish between expressions, variables, and numbers. Mathematically expressions have an "outer" function, whereas variables and numbers can be directly differentiated. The `isexpr` function in `TermInterface` returns `true` when passed an expression, and `false` when passed a symbol or numeric literal. The latter two may be distinguished by `isa(..., Symbol)`. @@ -75,8 +66,8 @@ Here we create a function, `D`, that when it encounters an expression it *dispat function D(ex, var=:x) if isexpr(ex) op, args = operation(ex), arguments(ex) - D(Val(op), arity(ex), args, var) - elseif isa(ex, Symbol) && ex == :x + D(Val(op), Val(arity(ex)), args, var) + elseif isa(ex, Symbol) && ex == var 1 else 0 @@ -84,6 +75,8 @@ function D(ex, var=:x) end ``` +(The use of `Val` is an idiom of `Julia` allowing dispatch on certain values such as function names and numbers.) + Now to develop methods for `D` for different "outside" functions and arities. @@ -91,31 +84,31 @@ Addition can be unary (`:(+x)` is a valid quoting, even if it might simplify to ```{julia} -D(::Val{:+}, ::Val{:unary}, args, var) = D(first(args), var) +D(::Val{:+}, ::Val{1}, args, var) = D(first(args), var) -function D(::Val{:+}, ::Val{:binary}, args, var) +function D(::Val{:+}, ::Val{2}, args, var) aβ€², bβ€² = D.(args, var) :($aβ€² + $bβ€²) end -function D(::Val{:+}, ::Val{:nary}, args, var) +function D(::Val{:+}, ::Any, args, var) aβ€²s = D.(args, var) :(+($aβ€²s...)) end ``` -The `args` are always held in a container, so the unary method must pull out the first one. The binary case should read as: apply `D` to each of the two arguments, and then create a quoted expression containing the sum of the results. The dollar signs interpolate into the quoting. (The "primes" are unicode notation achieved through `\prime[tab]` and not operations.) The *nary* case does something similar, only using splatting to produce the sum. +The `args` are always held in a container, so the unary method must pull out the first one. The binary case should read as: apply `D` to each of the two arguments, and then create a quoted expression containing the sum of the results. The dollar signs interpolate into the quoting. (The "primes" are unicode notation achieved through `\prime[tab]` and not operations.) The *nary* method (which catches *any* arity besides `1` and `2`) does something similar, only using splatting to produce the sum. -Subtraction must also be implemented in a similar manner, but not for the *nary* case: +Subtraction must also be implemented in a similar manner, but not for the *nary* case, as subtraction is not associative: ```{julia} -function D(::Val{:-}, ::Val{:unary}, args, var) +function D(::Val{:-}, ::Val{1}, args, var) aβ€² = D(first(args), var) :(-$aβ€²) end -function D(::Val{:-}, ::Val{:binary}, args, var) +function D(::Val{:-}, ::Val{2}, args, var) aβ€², bβ€² = D.(args, var) :($aβ€² - $bβ€²) end @@ -125,15 +118,15 @@ The *product rule* is similar to addition, in that $3$ cases are considered: ```{julia} -D(op::Val{:*}, ::Val{:unary}, args, var) = D(first(args), var) +D(op::Val{:*}, ::Val{1}, args, var) = D(first(args), var) -function D(::Val{:*}, ::Val{:binary}, args, var) +function D(::Val{:*}, ::Val{2}, args, var) a, b = args aβ€², bβ€² = D.(args, var) :($aβ€² * $b + $a * $bβ€²) end -function D(op::Val{:*}, ::Val{:nary}, args, var) +function D(op::Val{:*}, ::Any, args, var) a, bs... = args b = :(*($(bs...))) aβ€² = D(a, var) @@ -149,7 +142,7 @@ Division is only a binary operation, so here we have the *quotient rule*: ```{julia} -function D(::Val{:/}, ::Val{:binary}, args, var) +function D(::Val{:/}, ::Val{2}, args, var) u,v = args uβ€², vβ€² = D(u, var), D(v, var) :( ($uβ€²*$v - $u*$vβ€²)/$v^2 ) @@ -160,7 +153,7 @@ Powers are handled a bit differently. The power rule would require checking if t ```{julia} -function D(::Val{:^}, ::Val{:binary}, args, var) +function D(::Val{:^}, ::Val{2}, args, var) a, b = args D(:(exp($b*log($a))), var) # a > 0 assumed here end @@ -170,26 +163,26 @@ That leaves the task of defining a rule to differentiate both `exp` and `log`. W ```{julia} -function D(::Val{:exp}, ::Val{:unary}, args, var) - a = first(args) +function D(::Val{:exp}, ::Val{1}, args, var) + a = only(args) aβ€² = D(a, var) :(exp($a) * $aβ€²) end -function D(::Val{:log}, ::Val{:unary}, args, var) - a = first(args) +function D(::Val{:log}, ::Val{1}, args, var) + a = only(args) aβ€² = D(a, var) :(1/$a * $aβ€²) end -function D(::Val{:sin}, ::Val{:unary}, args, var) - a = first(args) +function D(::Val{:sin}, ::Val{1}, args, var) + a = only(args) aβ€² = D(a, var) :(cos($a) * $aβ€²) end -function D(::Val{:cos}, ::Val{:unary}, args, var) - a = first(args) +function D(::Val{:cos}, ::Val{1}, args, var) + a = only(args) aβ€² = D(a, var) :(-sin($a) * $aβ€²) end @@ -207,16 +200,16 @@ More functions could be included, but for this example the above will suffice, a ```{julia} -ex₁ = :(x + 2/x) -D(ex₁, :x) +ex = :(x + 2/x) +D(ex, :x) ``` The output does not simplify, so some work is needed to identify `1 - 2/x^2` as the answer. ```{julia} -exβ‚‚ = :( (x + sin(x))/sin(x)) -D(exβ‚‚, :x) +ex = :( (x + sin(x))/sin(x)) +D(ex, :x) ``` Again, simplification is not performed. @@ -226,8 +219,8 @@ Finally, we have a second derivative taken below: ```{julia} -ex₃ = :(sin(x) - x - x^3/6) -D(D(ex₃, :x), :x) +ex = :(sin(x) - x - x^3/6) +D(D(ex, :x), :x) ``` The length of the expression should lead to further appreciation for simplification steps taken when doing such a computation by hand. diff --git a/quarto/index.qmd b/quarto/index.qmd index 6453f08..8d3ab94 100644 --- a/quarto/index.qmd +++ b/quarto/index.qmd @@ -72,7 +72,7 @@ books visit . These notes may be compiled into a `pdf` file through Quarto. As the result is rather large, we do not provide that file for download. For the interested reader, downloading the repository, instantiating the environment, and running `quarto` to render to `pdf` in the `quarto` subdirectory should produce that file (after some time). --> -To *contribute* -- say by suggesting addition topics, correcting a +To *contribute* -- say by suggesting additional topics, correcting a mistake, or fixing a typo -- click the "Edit this page" link and join the list of [contributors](https://github.com/jverzani/CalculusWithJuliaNotes.jl/graphs/contributors). Thanks to all contributors and a *very* special thanks to `@fangliu-tju` for their careful and most-appreciated proofreading. ## Running Julia diff --git a/quarto/references.bib b/quarto/references.bib index 7c85495..48a8b35 100644 --- a/quarto/references.bib +++ b/quarto/references.bib @@ -63,3 +63,27 @@ pages = {97–111}, numpages = {15} } + +@book{doi:10.1137/1.9781611977165, +author = {Trefethen, Lloyd N. and Bau, David}, +title = {Numerical Linear Algebra, Twenty-fifth Anniversary Edition}, +publisher = {Society for Industrial and Applied Mathematics}, +year = {2022}, +doi = {10.1137/1.9781611977165}, +address = {Philadelphia, PA}, +edition = {}, +URL = {https://epubs.siam.org/doi/abs/10.1137/1.9781611977165}, +eprint = {https://epubs.siam.org/doi/pdf/10.1137/1.9781611977165} +} + +@book{doi:10.1137/1.9781611975086, +author = {Driscoll, Tobin A. and Braun, Richard J.}, +title = {Fundamentals of Numerical Computation}, +publisher = {Society for Industrial and Applied Mathematics}, +year = {2017}, +doi = {10.1137/1.9781611975086}, +address = {Philadelphia, PA}, +edition = {Julia adaptation}, +URL = {https://tobydriscoll.net/fnc-julia/frontmatter.html}, +eprint = {https://epubs.siam.org/doi/pdf/10.1137/1.9781611975086} +}