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

Merge branch 'main' into ranges

Browse Source
This commit is contained in:
Andrew Lin
2022-08-23 08:54:45 -07:00
parent 93761d8c24 c9f10fe2d7
commit 6517e7fb7d
12 changed files with 1467 additions and 424 deletions
Download Patch File Download Diff File Expand all files Collapse all files

392
CwJ/alternatives/symbolics.__jmd__
View File

@@ -1,392 +0,0 @@
# Symbolics.jl
Incorporate:
Basics
https://github.com/SciML/ModelingToolkit.jl
https://github.com/JuliaSymbolics/Symbolics.jl
https://github.com/JuliaSymbolics/SymbolicUtils.jl
* Rewriting
https://github.com/JuliaSymbolics/SymbolicUtils.jl
* Plotting
Polynomials
Limits
XXX ... room here!
Derivatives
https://github.com/JuliaSymbolics/Symbolics.jl
Integration
https://github.com/SciML/SymbolicNumericIntegration.jl
The `Symbolics.jl` package is a Computer Algebra System (CAS) built entirely in `Julia`.
This package is under heavy development.
## Algebraic manipulations
### construction
@variables
SymbolicUtils.@syms assumptions
x is a `Num`, `Symbolics.value(x)` is of type `SymbolicUtils{Real, Nothing}
relation to SymbolicUtils
Num wraps things; Term
### Substitute
### Simplify
simplify
expand
rewrite rules
### Solving equations
solve_for
## Expressions to functions
build_function
## Derivatives
1->1: Symbolics.derivative(x^2 + cos(x), x)
1->3: Symbolics.derivative.([x^2, x, cos(x)], x)
3 -> 1: Symbolics.gradient(x*y^z, [x,y,z])
2 -> 2: Symbolics.jacobian([x,y^z], [x,y])
# higher order
1 -> 1: D(ex, x, n=1) = foldl((ex,_) -> Symbolics.derivative(ex, x), 1:n, init=ex)
2 -> 1: (2nd) Hessian
## Differential equations
## Integrals
WIP
## ----
# follow sympy tutorial
using Symbolics
import SymbolicUtils
@variables x y z
# substitution
ex = cos(x) + 1
substitute(ex, Dict(x=>y))
substitute(ex, Dict(x=>0)) # does eval
ex = x^y
substitute(ex, Dict(y=> x^y))
# expand trig
r1 = @rule sin(2 * ~x) => 2sin(~x)*cos(~x)
r2 = @rule cos(2 * ~x) => cos(~x)^2 - sin(~x)^2
expand_trig(ex) = simplify(ex, RuleSet([r1, r2]))
ex = sin(2x) + cos(2x)
expand_trig(ex)
## Multiple
@variables x y z
ex = x^3 + 4x*y -z
substitute(ex, Dict(x=>2, y=>4, z=>0))
# Converting Strings to Expressions
# what is sympify?
# evalf
# lambdify: symbolic expression -> function
ex = x^3 + 4x*y -z
λ = build_function(ex, x,y,z, expression=Val(false))
λ(2,4,0)
# pretty printing
using Latexify
latexify(ex)
# Simplify
@variables x y z t
simplify(sin(x)^2 + cos(x)^2)
simplify((x^3 + x^2 - x - 1) / (x^2 + 2x + 1)) # fails, no factor
simplify(((x+1)*(x^2-1))/((x+1)^2)) # works
import SpecialFunctions: gamma
simplify(gamma(x) / gamma(x-2)) # fails
# Polynomial
## expand
expand((x+1)^2)
expand((x+2)*(x-3))
expand((x+1)*(x-2) - (x-1)*x)
## factor
### not defined
## collect
COLLECT_RULES = [
@rule(~x*x^(~n::SymbolicUtils.isnonnegint) => (~x, ~n))
@rule(~x * x => (~x, 1))
]
function _collect(ex, x)
d = Dict()
exs = expand(ex)
if SymbolicUtils.operation(Symbolics.value(ex)) != +
d[0] => ex
else
for aᵢ ∈ SymbolicUtils.arguments(Symbolics.value(expand(ex)))
u = simplify(aᵢ, RuleSet(COLLECT_RULES))
if isa(u, Tuple)
a,n = u
else
a,n = u,0
end
d[n] = get(d, n, 0) + a
end
end
d
end
## cancel -- no factor
## apart -- no factor
## Trignometric simplification
INVERSE_TRIG_RUELS = [@rule(cos(acos(~x)) => ~x)
@rule(acos(cos(~x)) => abs(rem2pi(~x, RoundNearest)))
@rule(sin(asin(~x)) => ~x)
@rule(asin(sin(~x)) => abs(rem2pi(x + pi/2, RoundNearest)) - pi/2)
]
@variables θ
simplify(cos(acos(θ)), RuleSet(INVERSE_TRIG_RUELS))
# Copy from https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/master/src/simplify_rules.jl
# the TRIG_RULES are applied by simplify by default
HTRIG_RULES = [
@acrule(-sinh(~x)^2 + cosh(~x)^2 => one(~x))
@acrule(sinh(~x)^2 + 1 => cosh(~x)^2)
@acrule(cosh(~x)^2 + -1 => -sinh(~x)^2)
@acrule(tanh(~x)^2 + 1*sech(~x)^2 => one(~x))
@acrule(-tanh(~x)^2 + 1 => sech(~x)^2)
@acrule(sech(~x)^2 + -1 => -tanh(~x)^2)
@acrule(coth(~x)^2 + -1*csch(~x)^2 => one(~x))
@acrule(coth(~x)^2 + -1 => csch(~x)^2)
@acrule(csch(~x)^2 + 1 => coth(~x)^2)
@acrule(tanh(~x) => sinh(~x)/cosh(~x))
@acrule(sinh(-~x) => -sinh(~x))
@acrule(cosh(-~x) => -cosh(~x))
]
trigsimp(ex) = simplify(simplify(ex, RuleSet(HTRIG_RULES)))
trigsimp(sin(x)^2 + cos(x)^2)
trigsimp(sin(x)^4 -2cos(x)^2*sin(x)^2 + cos(x)^4) # no factor
trigsimp(cosh(x)^2 + sinh(x)^2)
trigsimp(sinh(x)/tanh(x))
EXPAND_TRIG_RULES = [
@acrule(sin(~x+~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y))
@acrule(sinh(~x+~y) => sinh(~x)*cosh(~y) + cosh(~x)*sinh(~y))
@acrule(sin(2*~x) => 2sin(~x)*cos(~x))
@acrule(sinh(2*~x) => 2sinh(~x)*cosh(~x))
@acrule(cos(~x+~y) => cos(~x)*cos(~y) - sin(~x)*sin(~y))
@acrule(cosh(~x+~y) => cosh(~x)*cosh(~y) + sinh(~x)*sinh(~y))
@acrule(cos(2*~x) => cos(~x)^2 - sin(~x)^2)
@acrule(cosh(2*~x) => cosh(~x)^2 + sinh(~x)^2)
@acrule(tan(~x+~y) => (tan(~x) - tan(~y)) / (1 + tan(~x)*tan(~y)))
@acrule(tanh(~x+~y) => (tanh(~x) + tanh(~y)) / (1 + tanh(~x)*tanh(~y)))
@acrule(tan(2*~x) => 2*tan(~x)/(1 - tan(~x)^2))
@acrule(tanh(2*~x) => 2*tanh(~x)/(1 + tanh(~x)^2))
]
expandtrig(ex) = simplify(simplify(ex, RuleSet(EXPAND_TRIG_RULES)))
expandtrig(sin(x+y))
expandtrig(tan(2x))
# powers
# in genearl x^a*x^b = x^(a+b)
@variables x y a b
simplify(x^a*x^b - x^(a+b)) # 0
# x^a*y^a = (xy)^a When x,y >=0, a ∈ R
simplify(x^a*y^a - (x*y)^a)
## ??? How to specify such assumptions?
# (x^a)^b = x^(ab) only if b ∈ Int
@syms x a b
simplify((x^a)^b - x^(a*b))
@syms x a b::Int
simplify((x^a)^b - x^(a*b)) # nope
ispositive(x) = isa(x, Real) && x > 0
_isinteger(x) = isa(x, Integer)
_isinteger(x::SymbolicUtils.Sym{T,S}) where {T <: Integer, S} = true
POWSIMP_RULES = [
@acrule((~x::ispositive)^(~a::isreal) * (~y::ispositive)^(~a::isreal) => (~x*~y)^~a)
@rule(((~x)^(~a))^(~b::_isinteger) => ~x^(~a * ~b))
]
powsimp(ex) = simplify(simplify(ex, RuleSet(POWSIMP_RULES)))
@syms x a b::Int
simplify((x^a)^b - x^(a*b)) # nope
EXPAND_POWER_RULES = [
@rule((~x)^(~a + ~b) => (_~)^(~a) * (~x)^(~b))
@rule((~x*~y)^(~a) => (~x)^(~a) * (~y)^(~a))
## ... more on simplification...
## Calculus
@variables x y z
import Symbolics: derivative
derivative(cos(x), x)
derivative(exp(x^2), x)
# multiple derivative
Symbolics.derivative(ex, x, n::Int) = reduce((ex,_) -> derivative(ex, x), 1:n, init=ex) # helper
derivative(x^4, x, 3)
ex = exp(x*y*z)
using Chain
@chain ex begin
derivative(x, 3)
derivative(y, 3)
derivative(z, 3)
end
# using Differential operator
expr = exp(x*y*z)
expr |> Differential(x)^2 |> Differential(y)^3 |> expand_derivatives
# no integrate
# no limit
# Series
function series(ex, x, x0=0, n=5)
Σ = zero(ex)
for i ∈ 0:n
ex = expand_derivatives((Differential(x))(ex))
Σ += substitute(ex, Dict(x=>0)) * x^i / factorial(i)
end
Σ
end
# finite differences
# Solvers
@variables x y z a
eq = x ~ a
Symbolics.solve_for(eq, x)
eqs = [x + y + z ~ 1
x + y + 2z ~ 3
x + 2y + 3z ~ 3
]
vars = [x,y,z]
xs = Symbolics.solve_for(eqs, vars)
[reduce((ex, r)->substitute(ex, r), Pair.(vars, xs), init=ex.lhs) for ex ∈ eqs] == [eq.rhs for eq ∈ eqs]
A = [1 1; 1 2]
b = [1, 3]
xs = Symbolics.solve_for(A*[x,y] .~ b, [x,y])
A*xs - b
A = [1 1 1; 1 1 2]
b = [1,3]
A*[x,y,z] - b
Symbolics.solve_for(A*[x,y,z] .~ b, [x,y,z]) # fails, singular
# nonlinear solve
# use `λ = mk_function(ex, args, expression=Val(false))`
# polynomial roots
# differential equations

658
CwJ/alternatives/symbolics.jmd Normal file
View File

@@ -0,0 +1,658 @@
# Symbolics.jl
There are a few options in `Julia` for symbolic math, for example, the `SymPy` package which wraps a Python library. This section describes a collection of native `Julia` packages providing many features of symbolic math.
## About
The `Symbolics` package bill itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which we don't describe here.
We begin by loading the `Symbolics` package which when loaded re-exports the `SymbolicUtils` package.
```julia
using Symbolics
```
## Symbolic variables
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two package discussed here.
`SymbolicUtils` creates variables which carry `Julia` type information (e.g. `Int`, `Float64`, ...). This type information carries through operations involving these variables. Symbolic variables can be created with the `@syms` macro. For example
```julia
@syms x y::Int f(x::Real)::Real
```
This creates `x` a symbolic value with type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputing a real variable.
The `symtype` function reveals the underlying type:
```julia
import Symbolics.SymbolicUtils: symtype
symtype(x), symtype(y)
```
For `y`, the symbolic type being real does not imply the `y` has a subtype of `Real`:
```julia
isa(y, Real)
```
We see that the function `f`, when called with `y` would return a value of (symbolic) type `Real`:
```julia
f(y) |> symtype
```
As the symbolic type of `x` is `Number` -- which is not a subtype of `Real` -- the following will error:
```julia; error=true
f(x)
```
!!! note
The `SymPy` package also has an `@syms` macro to create variables. Though they names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
For `Symbolics`, symbolic variables are created using a wrapper around an underlying `SymbolicUtils` object. This wrapper, `Num`, is a subtype of `Real` (the underlying `SymbolicUtils` object may have symbolic type `Real`, but it won't be a subtype of `Real`.)
Symbolic values are created with the `@variables` macro. For example:
```julia
@variables x y::Int z[1:3]::Int f(..)::Int
```
This creates
* a symbolic value `x` of `symtype` `Real`
* a symbolic value `y` of `symtype` `Int`
* a vector of symbolic values each of `symtype` `Int`
* a symbolic function `f` returning an object of `symtype` `Int`
The symbolic type reflects that of the underlying object behind the `Num` wrapper:
```julia
typeof(x), symtype(x), typeof(Symbolics.value(x))
```
(The `value` method unwraps the `Num` wrapper.)
## Symbolic expressions
Symbolic expressions are built up from symbolic variables through natural `Julia` idioms. `SymbolicUtils` privileges a few key operations: `Add`, `Mul`, `Pow`, and `Div`. For examples:
```julia
@syms x y
typeof(x + y) # `Add`
```
```julia
typeof(x * y) # `Mul`
```
Whereas, applying a function leaves a different type:
```julia
typeof(sin(x))
```
The `Term` wrapper just represents the effect of calling a function (in this case `sin`) on its arguments (in this case `x`).
This happens in the background with symbolic variables in `Symbolics`:
```julia
@variables x
typeof(sin(x)), typeof(Symbolics.value(sin(x)))
```
### Tree structure to expressions
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structdure of an expression. The main methods are (cf. [symbolicutils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
* `istree(ex)`: `true` if `ex` is not a *leaf* node (like a symbol or numeric literal)
* `operation(ex)`: the function being called (if `istree` returns `true`)
* `arguments(ex)`: the arguments to the function begin called
* `symtype(ex)`: the inferred type of the expression
In addition, the `issym` function, to determine if `x` is of type `Sym`, is useful to distinguish *leaf* nodes, as will be illustrated below.
These methods can be used to "walk" the tree:
```julia
@syms x y
ex = 1 + x^2 + y
operation(ex) # the outer function is `+`
```
```julia
arguments(ex) # `+` is n-ary, in this case with 3 arguments
```
```julia
ex1 = arguments(ex)[3] # terms have been reordered
operation(ex1) # operation for `x^2` is `^`
```
```julia
a, b = arguments(ex1)
```
```julia
istree(ex1), istree(a)
```
Here `a` is not a "tree", as it has no operation or arguments, it is just a variable (the `x` variable).
The value of `symtype` is the *inferred* type of an expression, which may not match the actual type. For example,
```julia
@variables x::Int
symtype(x), symtype(sin(x)), symtype(x/x), symtype(x / x^2)
```
The last one, is not likely to be an integer, but that is the inferred type in this case.
##### Example
As an example, we write a function to find the free symbols in a symbolic expression comprised of `SymbolicUtils` variables. (The `Symbolics.get_variables` also does this task.) To find the symbols involves walking the expression tree until a leaf node is found and then adding that to our collection if it matches `issym`.
```julia
import Symbolics.SymbolicUtils: issym
free_symbols(ex) = (s=Set(); free_symbols!(s, ex); s)
function free_symbols!(s, ex)
if istree(ex)
for a ∈ arguments(ex)
free_symbols!(s, a)
end
else
issym(ex) && push!(s, ex) # push new symbol onto set
end
end
```
```julia
@syms x y z
ex = sin(x + 1)*cos(z)
free_symbols(ex)
```
## Expression manipulation
### Substitute
The `substitute` command is used to replace values with other values. For examples:
```julia
@variables x y z
ex = 1 + x + x^2/2 + x^3/6
substitute(ex, x=>1)
```
This defines a symbolic expression, then substitutes the value `1` in for `x`. The `Pair` notation is useful for a *single* substitution. When there is more than one substitution, a dictionary is used:
```julia
w = x^3 + y^3 - 2z^3
substitute(w, Dict(x=>2, y=>3))
```
The `fold` argument can be passed `false` to inhibit evaluation of values. Compare:
```julia
ex = 1 + sqrt(x)
substitute(ex, x=>2), substitute(ex, x=>2, fold=false)
```
Or
```julia
ex = sin(x)
substitute(ex, x=>π), substitute(ex, x=>π, fold=false)
```
### Simplify
Algebraic operations with symbolic values can involve an exponentially increasing number of terms. As such, some simplification rules are applied after an operation to reduce the complexity of the computed value.
For example, `0+x` should simplify to `x`, `0+x`, `x^0`, or `x^1` should simplify, to some natural answer.
`SymbolicUtils` also [simplifies](https://symbolicutils.juliasymbolics.org/#simplification) several other expressions, including:
* `-x` becomes `(-1)*x`
* `x * x` becomes `x^2` (and `x^n` if more terms). Meaning this expression is represented as a power, not a product
* `x + x` becomes `2*x` (and `n*x` if more terms). Similarly, this represented as a product, not a sum.
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`
In `SymbolicUtils`, this *rewriting* is accomplished by means of *rewrite rules*. The package makes it easy to apply user-written rewrite rules.
### Rewriting
Many algebraic simplifications are done by the `simplify` command. For example, the basic trignometric identities are applied:
```julia
@variables x
ex = sin(x)^2 + cos(x)^2
ex, simplify(ex)
```
The `simplify` function applies a series of rewriting rule until the expression stabilizes. The rewrite rules can be user generated, if desired. For example, the Pythagorean identity of trigonometry, just used, can be implement with this rule:
```julia
r = @acrule(sin(~x)^2 + cos(~x)^2 => one(~x))
ex |> Symbolics.value |> r |> Num
```
The rewrite rule, `r`, is defined by the `@acrule` macro. The `a` is for associative, the `c` for commutative, assumptions made by the macro. (The `c` means `cos(x)^2 + sin(x)^2` will also simplify.) Rewrite rules are called on the underlying `SymbolicUtils` expression, so we first unwrap, then after re-wrap.
The above expression for `r` is fairly easy to appreciate. The value `~x` matches the same variable or expression. So the above rule will also simplify more complicated expressions:
```julia
@variables y z
ex1 = substitute(ex, x => sin(x + y + z))
ex1 |> Symbolics.value |> r |> Num
```
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trignometry and exponential simplifications:
```julia
r = @rule(exp(~x)^(~y) => exp(~x * ~y))
ex = exp(-x+z)^y
ex, ex |> Symbolics.value |> r |> Num
```
This rule is not commutative or associative, as `x^y` is not the same as `y^x` and `(x^y)^z` is not `x^(y^z)` in general.
The application of rules can be filtered through qualifying predicates. This artificial example uses `iseven` which returns `true` for even numbers. Here we subtract `1` when a number is not even, and otherwise leave the number alone. We do this with two rules:
```julia
reven = @rule ~x::iseven => ~x
rodd = @rule ~x::(!iseven) => ~x - 1
r = SymbolicUtils.Chain([rodd, reven])
r(2), r(3)
```
The `Chain` function conveniently allows the sequential application of rewrite rules.
The notation `~x` is called a "slot variable" in the [documentation](https://symbolicutils.juliasymbolics.org/rewrite/) for `SymbolicUtils`. It matches a single expression. To match more than one expression, a "segment variable", denoted with two `~`s is used.
### Creating functions
By utilizing the tree-like nature of a symbolic expression, a `Julia` expression can be built from an symbolic expression easily enough. The `Symbolics.toexpr` function does this:
```julia
ex = exp(-x + z)^y
Symbolics.toexpr(ex)
```
This output shows an internal representation of the steps for computing the value `ex` given different inputs.
Such `Julia` expressions are one step away from building `Julia` functions for evaluating symbolic expressions fast (though with some technical details about "world age" to be reckoned with). The `build_function` function with the argument `expression=Val(false)` will compile a `Julia` function:
```julia
h = build_function(ex, x, y, z; expression=Val(false))
h(1, 2, 3)
```
The above is *similar* to substitution:
```julia
substitute(ex, Dict(x=>1, y=>2, z=>3))
```
However, `build_function` will be **significantly** more performant, which when many function calls are used -- such as with plotting -- is a big advantage.
!!! note
The documentation colorfully says "`build_function` is kind of like if `lambdify` (from `SymPy`) ate its spinach."
The above, through passing ``3`` variables after the expression, creates a function of ``3`` variables. Functions of a vector of inputs can also be created, just by expressing the variables in that manner:
```juila
h1 = build_function(ex, [x, y, z]; expression=Val(false))
h1([1, 2, 3]) # not h1(1,2,3)
```
##### Example
As an example, here we use the `Roots` package to find a zero of a function defined symbolically:
```julia
import Roots
@variables x
ex = x^5 - x - 1
λ = build_function(ex, x; expression=Val(false))
Roots.find_zero(λ, (1, 2))
```
### Plotting
Using `Plots`, the plotting symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
For example,
```julia
using Plots
@variables x
plot(x^x^x, 0, 2)
```
A parametric plot is easily defined:
```julia
plot(sin(x), cos(x), 0, pi/4)
```
Expressions to be plotted can represent multivariate functions.
```julia
@variables x y
ex = 3*(1-x)^2*exp(-x^2 - (y+1)^2) - 10(x/5-x^3-y^5)*exp(-x^2-y^2) - 1/3*exp(-(x+1)^2-y^2)
xs = ys = range(-5, 5, length=100)
surface(xs, ys, ex)
```
The ordering of the variables is determined by `Symbolics.get_variables`:
```julia
Symbolics.get_variables(ex)
```
### Polynomial manipulations
There are some facilities for manipulating polynomial expressions in `Symbolics`. A polynomial, mathematically, is an expression involving one or more symbols with coefficients from a collection that has, at a minimum, addition and multiplication defined. The basic building blocks of polynomials are *monomials*, which are comprised of products of powers of the symbols. Mathematically, monomials are often allowed to have a multiplying coefficient and may be just a coefficient (if each symbol is taken to the power ``0``), but here we consider just expressions of the type ``x_1^{a_1} \cdot x_2^{a_2} \cdots x_k^{a_k}`` with the ``a_i > 0`` as monomials.
With this understanding, then an expression can be broken up into monomials with a possible leading coefficient (possibly ``1``) *and* terms which are not monomials (such as a constant or a more complicated function of the symbols). This is what is returned by the `polynomial_coeffs` function.
For example
```julia
@variables a b c x
d, r = polynomial_coeffs(a*x^2 + b*x + c, (x,))
```
The first term output is dictionary who's keys are the monomials and who's values are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The expression can then be reconstructed through
```julia
r + sum(v*k for (k,v) ∈ d)
```
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomials symbols to `polynomial_coeff` as an iterable. (Above as a ``1``-element tuple.)
More complicated polynomials can be similarly decomposed:
```julia
@variables a b c x y z
ex = a*x^2*y*z + b*x*y^2*z + c*x*y*z^2
d, r = polynomial_coeffs(ex, (x, y, z))
```
The (sparse) decomposition of the polynomial is returned through `d`. The same pattern as above can be used to reconstruct the expression.
To extract the coefficient for a monomial term, indexing can be used. Of note, is an expression like `x^2*y*z` could *possibly* not equal the algebraically equal `x*y*z*x`, as they are only equal after some simplification, but the keys are in a canonical form, so this is not a concern:
```julia
d[x*y*z*x], d[z*y*x^2]
```
The residual term will capture any non-polynomial terms:
```julia
ex = sin(x) - x + x^3/6
d, r = polynomial_coeffs(ex, (x,))
r
```
To find the degree of a monomial expression, the `degree` function is available. Here it is applied to each monomial in `d`:
```julia
[degree(k) for (k,v) ∈ d]
```
The `degree` function will also identify the degree of more complicated terms:
```julia
degree(1 + x + x^2)
```
Mathematically the degree of the ``0`` polynomial may be ``-1`` or undefined, but here it is ``0``:
```julia
degree(0), degree(1), degree(x), degree(x^a)
```
The coefficients are returned as *values* of a dictionary, and dictionaries are unsorted. To have a natural map between polynomials of a single symbol in the standard basis and a vector, we could use a pattern like this:
```julia
@variables x a0 as[1:10]
p = a0 + sum(as[i]*x^i for i ∈ eachindex(collect(as)))
d, r = polynomial_coeffs(p, (x,))
d
```
To sort the values we can use a pattern like the following:
```julia
vcat(r, [d[k] for k ∈ sort(collect(keys(d)), by=degree)])
```
----
Rational expressions can be decomposed into a numerator and denominator using the following idiom, which ensures the outer operation is division (a binary operation):
```julia
@variables x
ex = (1 + x + x^2) / (1 + x + x^2 + x^3)
function nd(ex)
ex1 = Symbolics.value(ex)
(operation(ex1) == /) || return (ex, one(ex))
Num.(arguments(ex1))
end
nd(ex)
```
With this, the study of asymptotic behaviour of a univariate rational expression would involve an investigation like the following:
```julia
m,n = degree.(nd(ex))
m > n ? "limit is infinite" : m < n ? "limit is 0" : "limit is a constant"
```
### Algebraically solving equations
The `~` operator creates a symbolic equation. For example
```julia
@variables x y
x^5 - x ~ 1
```
or
```julia
ex = [5x + 2y, 6x + 3y] .~ [1, 2]
```
The `Symbolics.solve_for` function can solve *linear* equations. For example,
```julia
Symbolics.solve_for(ex, [x, y])
```
### Limits
As of writing, there is no extra functionality provided by `Symbolics` for computing limits.
### Derivatives
`Symbolics` provides the `derivative` function to compute the derivative of a function with respect to a variable:
```julia
@variables a b c x
ex = a*x^2 + b*x + c
Symbolics.derivative(ex, x)
```
The computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
```julia
D = Differential(x)
D(ex)
```
and then
```julia
expand_derivatives(D(ex))
```
The differentials can be multiplied to create operators for taking higher-order derivatives:
```julia
@variables x y
ex = (x - y^2)/(x^2 + y^2)
Dx, Dy = Differential(x), Differential(y)
Dxx, Dxy, Dyy = Dx*Dx, Dx*Dy, Dy*Dy
[Dxx(ex) Dxy(ex); Dxy(ex) Dyy(ex)] .|> expand_derivatives
```
In addition to `Symbolics.derivative` there are also the helper functions, such as `hessian` which performs the above
```julia
Symbolics.hessian(ex, [x,y])
```
The `gradient` function is also available
```julia
@variables x y z
ex = x^2 - 2x*y + z*y
Symbolics.gradient(ex, [x, y, z])
```
The `jacobian` takes an array of expressions:
```julia
@variables x y
exs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(exs, [x,y])
```
### Integration
The `SymbolicNumericIntegration` package provides a means to integrate *univariate* expressions through its `integrate` function.
Symbolic integration can be approached in different ways. SymPy implements part of the Risch algorithm in addition to other algorithms. Rules-based algorithms could also be implemented.
For example, here is a simple rule that could be used to integrate a single integral
```julia
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
@syms x ∫(x)
r = @rule ∫(~x::is_var) => x^2/2
r(∫(x))
```
The `SymbolicNumericIntegration` package includes many more predicates for doing rules-based integration, but it primarily approaches the task in a different manner.
If ``f(x)`` is to be integrated, a set of *candidate* answers is generated. The following is **proposed** as an answer: ``\sum q_i \Theta_i(x)``. Differentiating the proposed answer leads to a *linear system of equations* that can be solved.
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with ``f(x) = x \sin(x)`` and the candidate thetas are ``{x, \sin(x), \cos(x), x\sin(x), x\cos(x)}`` so that we propose:
```math
\int f(x) dx = q_1 x + q_2 \sin(x) + q_3 \cos(x) + q_4 x \sin(x) + q_4 x \cos(x)
```
Differentiating both sides, yields a term ``x\sin(x)`` on the left, and equating coefficients gives:
```math
q_1 = q_4 = 0,\quad q_5 = -1, \quad q_4 - q_3 = q_2 - q_5 = 0
```
which can be solved with ``q_5=-1``, ``q_2=1``, and the other coefficients being ``0``. That is ``\int f(x) dx = 1 \sin(x) + (-1) x\cos(x)``.
The package provides an algorithm for the creation of candidates and the means to solve when possible. The `integrate` function is the main entry point. It returns three values: `solved`, `unsolved`, and `err`. The `unsolved` is the part of the integrand which can not be solved through this package. It is `0` for a given problem when `integrate` is successful in identifying an antiderivative, in which case `solved` is the answer. The value of `err` is a bound on the numerical error introduced by the algorithm.
To see, we have:
```julia
using SymbolicNumericIntegration
@variables x
integrate(x * sin(x))
```
The second term is `0`, as this has an identified antiderivative.
```julia
integrate(exp(x^2) + sin(x))
```
This returns `exp(x^2)` for the unsolved part, as this function has no simple antiderivative.
Powers of trig functions have antiderivatives, as can be deduced using integration by parts. When the fifth power is used, there is a numeric aspect to the algorithm that is seen:
```julia
u,v,w = integrate(sin(x)^5)
```
The derivative of `u` matches up to some numeric tolerance:
```julia
Symbolics.derivative(u, x) - sin(x)^5
```
The integration of rational functions (ratios of polynomials) can be done algorithmically, provided the underlying factorizations can be identified. The `SymbolicNumericIntegration` package has a function `factor_rational` that can identify factorizations.
```julia
import SymbolicNumericIntegration: factor_rational
@variables x
u = (1 + x + x^2)/ (x^2 -2x + 1)
v = factor_rational(u)
```
The summands in `v` are each integrable. We can see that `v` is a reexpression through
```julia
simplify(u - v)
```
The algorithm is numeric, not symbolic. This can be seen in these two factorizations:
```julia
u = 1 / expand((x^2-1)*(x-2)^2)
v = factor_rational(u)
```
or
```julia
u = 1 / expand((x^2+1)*(x-2)^2)
v = factor_rational(u)
```
As such, the integrals have numeric differences:
```julia
a,b,c = integrate(u)
```
We can see a bit of why through the following which needs a tolerance set to identify the rational numbers correctly:
```julia
cs = [first(arguments(term)) for term ∈ arguments(a)] # pick off coefficients
```
```julia
rationalize.(cs; tol=1e-8)
```

2
CwJ/integrals/ftc.jmd
View File

@@ -413,7 +413,7 @@ So the answer to the question is
This seems like a lot of work, and indeed it is more than is needed. The following would be more typical once the rules are learned:
```math
\int_0^\pi 100 \sin(x) dx = -100(-\cos(x)) \big|_0^{\pi} = 100 \cos(x) \big|_{\pi}^0 = 100(1) - 100(-1) = 200.
\int_0^\pi 100 \sin(x) dx = 100(-\cos(x)) \big|_0^{\pi} = 100 \cos(x) \big|_{\pi}^0 = 100(1) - 100(-1) = 200.
```
## The derivative of the integral

2
CwJ/precalc/polynomial.jmd
View File

@@ -180,7 +180,7 @@ manner. This allows the `Julia` package `SymPy` to provide
functionality from SymPy within `Julia`.
!!! note
When `SymPy` is installed through the package manger, the underlying `Python`
When `SymPy` is installed through the package manager, the underlying `Python`
libraries will also be installed.
!!! note

6
quarto/README-quarto.md
View File

@@ -33,7 +33,7 @@ To compile the pages through quarto
* This error
> fatal: 'gh-pages' is already checked out at '/Users/verzani/julia/CalculusWithJuliaNotes/quarto/f5611730'
was solved with
was solved with (https://waylonwalker.com/til/git-checkout-worktree/)
> git worktree remove f5611730
@@ -45,9 +45,9 @@ Eventually, if this workflow seems to be settled:
* deprecate .jmd files
* deprecate need to make "pluto friendly"
* do something with JSXGraph
DONE? * do something with JSXGraph
* figure out why PlotlyLight doesn't work
* move to not use CalculusWithJulia.WeaveSupport
* use an include file not the "hack" in jmd2qmd
DONE * use an include file not the "hack" in jmd2qmd
* modify sympy's show method
* take advantage of mermaid, ojs, bibliography, ...

6
quarto/_quarto.yml
View File

@@ -1,4 +1,4 @@
version: 0.3
version: 0.4
project:
type: book
@@ -108,9 +108,9 @@ book:
- integral_vector_calculus/stokes_theorem.qmd
- integral_vector_calculus/review.qmd
- part: "Alternatives"
- part: "Alternative packages"
chapters:
# - alternatives/symbolics.qmd
- alternatives/symbolics.qmd
# - alternatives/sciML.qmd
# - alternatives/interval_arithmetic.qmd
- alternatives/plotly_plotting.qmd

778
quarto/alternatives/symbolics.qmd Normal file
View File

@@ -0,0 +1,778 @@
# Symbolics.jl
There are a few options in `Julia` for symbolic math, for example, the `SymPy` package which wraps a Python library. This section describes a collection of native `Julia` packages providing many features of symbolic math.
## About
The `Symbolics` package bill itself as a "fast and modern Computer Algebra System (CAS) for a fast and modern programming language." This package relies on the `SymbolicUtils` package and is built upon by the `ModelingToolkit` package, which we don't describe here.
We begin by loading the `Symbolics` package which when loaded re-exports the `SymbolicUtils` package.
```{julia}
using Symbolics
```
## Symbolic variables
Symbolic math at its core involves symbolic variables, which essentially defer evaluation until requested. The creation of symbolic variables differs between the two package discussed here.
`SymbolicUtils` creates variables which carry `Julia` type information (e.g. `Int`, `Float64`, ...). This type information carries through operations involving these variables. Symbolic variables can be created with the `@syms` macro. For example
```{julia}
@syms x y::Int f(x::Real)::Real
```
This creates `x` a symbolic value with type `Number`, `y` a symbolic variable holding integer values, and `f` a symbolic function of a single real variable outputing a real variable.
The `symtype` function reveals the underlying type:
```{julia}
import Symbolics.SymbolicUtils: symtype
symtype(x), symtype(y)
```
For `y`, the symbolic type being real does not imply the `y` has a subtype of `Real`:
```{julia}
isa(y, Real)
```
We see that the function `f`, when called with `y` would return a value of (symbolic) type `Real`:
```{julia}
f(y) |> symtype
```
As the symbolic type of `x` is `Number` – which is not a subtype of `Real` – the following will error:
```{julia}
#| error: true
f(x)
```
:::{.callout-note}
## Note
The `SymPy` package also has an `@syms` macro to create variables. Though they names agree, they do different things. Using both packages together would require qualifying many shared method names. For `SymbolicUtils`, the `@syms` macro uses `Julia` types to parameterize the variables. In `SymPy` it is possible to specify *assumptions* on the variables, but that is different and not useful for dispatch without some extra effort.
:::
For `Symbolics`, symbolic variables are created using a wrapper around an underlying `SymbolicUtils` object. This wrapper, `Num`, is a subtype of `Real` (the underlying `SymbolicUtils` object may have symbolic type `Real`, but it won't be a subtype of `Real`.)
Symbolic values are created with the `@variables` macro. For example:
```{julia}
@variables x y::Int z[1:3]::Int f(..)::Int
```
This creates
* a symbolic value `x` of `symtype` `Real`
* a symbolic value `y` of `symtype` `Int`
* a vector of symbolic values each of `symtype` `Int`
* a symbolic function `f` returning an object of `symtype` `Int`
The symbolic type reflects that of the underlying object behind the `Num` wrapper:
```{julia}
typeof(x), symtype(x), typeof(Symbolics.value(x))
```
(The `value` method unwraps the `Num` wrapper.)
## Symbolic expressions
Symbolic expressions are built up from symbolic variables through natural `Julia` idioms. `SymbolicUtils` privileges a few key operations: `Add`, `Mul`, `Pow`, and `Div`. For examples:
```{julia}
@syms x y
typeof(x + y) # `Add`
```
```{julia}
typeof(x * y) # `Mul`
```
Whereas, applying a function leaves a different type:
```{julia}
typeof(sin(x))
```
The `Term` wrapper just represents the effect of calling a function (in this case `sin`) on its arguments (in this case `x`).
This happens in the background with symbolic variables in `Symbolics`:
```{julia}
@variables x
typeof(sin(x)), typeof(Symbolics.value(sin(x)))
```
### Tree structure to expressions
The `TermInterface` package is used by `SymbolicUtils` to explore the tree structdure of an expression. The main methods are (cf. [symbolicutils.jl](https://symbolicutils.juliasymbolics.org/#expression_interface)):
* `istree(ex)`: `true` if `ex` is not a *leaf* node (like a symbol or numeric literal)
* `operation(ex)`: the function being called (if `istree` returns `true`)
* `arguments(ex)`: the arguments to the function begin called
* `symtype(ex)`: the inferred type of the expression
In addition, the `issym` function, to determine if `x` is of type `Sym`, is useful to distinguish *leaf* nodes, as will be illustrated below.
These methods can be used to "walk" the tree:
```{julia}
@syms x y
ex = 1 + x^2 + y
operation(ex) # the outer function is `+`
```
```{julia}
arguments(ex) # `+` is n-ary, in this case with 3 arguments
```
```{julia}
ex1 = arguments(ex)[3] # terms have been reordered
operation(ex1) # operation for `x^2` is `^`
```
```{julia}
a, b = arguments(ex1)
```
```{julia}
istree(ex1), istree(a)
```
Here `a` is not a "tree", as it has no operation or arguments, it is just a variable (the `x` variable).
The value of `symtype` is the *inferred* type of an expression, which may not match the actual type. For example,
```{julia}
@variables x::Int
symtype(x), symtype(sin(x)), symtype(x/x), symtype(x / x^2)
```
The last one, is not likely to be an integer, but that is the inferred type in this case.
##### Example
As an example, we write a function to find the free symbols in a symbolic expression comprised of `SymbolicUtils` variables. (The `Symbolics.get_variables` also does this task.) To find the symbols involves walking the expression tree until a leaf node is found and then adding that to our collection if it matches `issym`.
```{julia}
import Symbolics.SymbolicUtils: issym
free_symbols(ex) = (s=Set(); free_symbols!(s, ex); s)
function free_symbols!(s, ex)
if istree(ex)
for a ∈ arguments(ex)
free_symbols!(s, a)
end
else
issym(ex) && push!(s, ex) # push new symbol onto set
end
end
```
```{julia}
@syms x y z
ex = sin(x + 1)*cos(z)
free_symbols(ex)
```
## Expression manipulation
### Substitute
The `substitute` command is used to replace values with other values. For examples:
```{julia}
@variables x y z
ex = 1 + x + x^2/2 + x^3/6
substitute(ex, x=>1)
```
This defines a symbolic expression, then substitutes the value `1` in for `x`. The `Pair` notation is useful for a *single* substitution. When there is more than one substitution, a dictionary is used:
```{julia}
w = x^3 + y^3 - 2z^3
substitute(w, Dict(x=>2, y=>3))
```
The `fold` argument can be passed `false` to inhibit evaluation of values. Compare:
```{julia}
ex = 1 + sqrt(x)
substitute(ex, x=>2), substitute(ex, x=>2, fold=false)
```
Or
```{julia}
ex = sin(x)
substitute(ex, x=>π), substitute(ex, x=>π, fold=false)
```
### Simplify
Algebraic operations with symbolic values can involve an exponentially increasing number of terms. As such, some simplification rules are applied after an operation to reduce the complexity of the computed value.
For example, `0+x` should simplify to `x`, `0+x`, `x^0`, or `x^1` should simplify, to some natural answer.
`SymbolicUtils` also [simplifies](https://symbolicutils.juliasymbolics.org/#simplification) several other expressions, including:
* `-x` becomes `(-1)*x`
* `x * x` becomes `x^2` (and `x^n` if more terms). Meaning this expression is represented as a power, not a product
* `x + x` becomes `2*x` (and `n*x` if more terms). Similarly, this represented as a product, not a sum.
* `p/q * x` becomes `(p*x)/q)`, similarly `p/q * x/y` becomes `(p*x)/(q*y)`
In `SymbolicUtils`, this *rewriting* is accomplished by means of *rewrite rules*. The package makes it easy to apply user-written rewrite rules.
### Rewriting
Many algebraic simplifications are done by the `simplify` command. For example, the basic trignometric identities are applied:
```{julia}
@variables x
ex = sin(x)^2 + cos(x)^2
ex, simplify(ex)
```
The `simplify` function applies a series of rewriting rule until the expression stabilizes. The rewrite rules can be user generated, if desired. For example, the Pythagorean identity of trigonometry, just used, can be implement with this rule:
```{julia}
r = @acrule(sin(~x)^2 + cos(~x)^2 => one(~x))
ex |> Symbolics.value |> r |> Num
```
The rewrite rule, `r`, is defined by the `@acrule` macro. The `a` is for associative, the `c` for commutative, assumptions made by the macro. (The `c` means `cos(x)^2 + sin(x)^2` will also simplify.) Rewrite rules are called on the underlying `SymbolicUtils` expression, so we first unwrap, then after re-wrap.
The above expression for `r` is fairly easy to appreciate. The value `~x` matches the same variable or expression. So the above rule will also simplify more complicated expressions:
```{julia}
@variables y z
ex1 = substitute(ex, x => sin(x + y + z))
ex1 |> Symbolics.value |> r |> Num
```
Rules involving two values are also easily created. This one, again, comes from the set of simplifications defined for trignometry and exponential simplifications:
```{julia}
r = @rule(exp(~x)^(~y) => exp(~x * ~y))
ex = exp(-x+z)^y
ex, ex |> Symbolics.value |> r |> Num
```
This rule is not commutative or associative, as `x^y` is not the same as `y^x` and `(x^y)^z` is not `x^(y^z)` in general.
The application of rules can be filtered through qualifying predicates. This artificial example uses `iseven` which returns `true` for even numbers. Here we subtract `1` when a number is not even, and otherwise leave the number alone. We do this with two rules:
```{julia}
reven = @rule ~x::iseven => ~x
rodd = @rule ~x::(!iseven) => ~x - 1
r = SymbolicUtils.Chain([rodd, reven])
r(2), r(3)
```
The `Chain` function conveniently allows the sequential application of rewrite rules.
The notation `~x` is called a "slot variable" in the [documentation](https://symbolicutils.juliasymbolics.org/rewrite/) for `SymbolicUtils`. It matches a single expression. To match more than one expression, a "segment variable", denoted with two `~`s is used.
### Creating functions
By utilizing the tree-like nature of a symbolic expression, a `Julia` expression can be built from an symbolic expression easily enough. The `Symbolics.toexpr` function does this:
```{julia}
ex = exp(-x + z)^y
Symbolics.toexpr(ex)
```
This output shows an internal representation of the steps for computing the value `ex` given different inputs.
Such `Julia` expressions are one step away from building `Julia` functions for evaluating symbolic expressions fast (though with some technical details about "world age" to be reckoned with). The `build_function` function with the argument `expression=Val(false)` will compile a `Julia` function:
```{julia}
h = build_function(ex, x, y, z; expression=Val(false))
h(1, 2, 3)
```
The above is *similar* to substitution:
```{julia}
substitute(ex, Dict(x=>1, y=>2, z=>3))
```
However, `build_function` will be **significantly** more performant, which when many function calls are used – such as with plotting – is a big advantage.
:::{.callout-note}
## Note
The documentation colorfully says "`build_function` is kind of like if `lambdify` (from `SymPy`) ate its spinach."
:::
The above, through passing $3$ variables after the expression, creates a function of $3$ variables. Functions of a vector of inputs can also be created, just by expressing the variables in that manner:
```{juila}
h1 = build_function(ex, [x, y, z]; expression=Val(false))
h1([1, 2, 3]) # not h1(1,2,3)
```
##### Example
As an example, here we use the `Roots` package to find a zero of a function defined symbolically:
```{julia}
import Roots
@variables x
ex = x^5 - x - 1
λ = build_function(ex, x; expression=Val(false))
Roots.find_zero(λ, (1, 2))
```
### Plotting
Using `Plots`, the plotting symbolic expressions is similar to the plotting of a function, as there is a plot recipe that converts the expression into a function via `build_function`.
For example,
```{julia}
using Plots
@variables x
plot(x^x^x, 0, 2)
```
A parametric plot is easily defined:
```{julia}
plot(sin(x), cos(x), 0, pi/4)
```
Expressions to be plotted can represent multivariate functions.
```{julia}
@variables x y
ex = 3*(1-x)^2*exp(-x^2 - (y+1)^2) - 10(x/5-x^3-y^5)*exp(-x^2-y^2) - 1/3*exp(-(x+1)^2-y^2)
xs = ys = range(-5, 5, length=100)
surface(xs, ys, ex)
```
The ordering of the variables is determined by `Symbolics.get_variables`:
```{julia}
Symbolics.get_variables(ex)
```
### Polynomial manipulations
There are some facilities for manipulating polynomial expressions in `Symbolics`. A polynomial, mathematically, is an expression involving one or more symbols with coefficients from a collection that has, at a minimum, addition and multiplication defined. The basic building blocks of polynomials are *monomials*, which are comprised of products of powers of the symbols. Mathematically, monomials are often allowed to have a multiplying coefficient and may be just a coefficient (if each symbol is taken to the power $0$), but here we consider just expressions of the type $x_1^{a_1} \cdot x_2^{a_2} \cdots x_k^{a_k}$ with the $a_i > 0$ as monomials.
With this understanding, then an expression can be broken up into monomials with a possible leading coefficient (possibly $1$) *and* terms which are not monomials (such as a constant or a more complicated function of the symbols). This is what is returned by the `polynomial_coeffs` function.
For example
```{julia}
@variables a b c x
d, r = polynomial_coeffs(a*x^2 + b*x + c, (x,))
```
The first term output is dictionary who's keys are the monomials and who's values are the coefficients. The second term, the residual, is all the remaining parts of the expression, in this case just the constant `c`.
The expression can then be reconstructed through
```{julia}
r + sum(v*k for (k,v) ∈ d)
```
The above has `a,b,c` as parameters and `x` as the symbol. This separation is designated by passing the desired polynomials symbols to `polynomial_coeff` as an iterable. (Above as a $1$-element tuple.)
More complicated polynomials can be similarly decomposed:
```{julia}
@variables a b c x y z
ex = a*x^2*y*z + b*x*y^2*z + c*x*y*z^2
d, r = polynomial_coeffs(ex, (x, y, z))
```
The (sparse) decomposition of the polynomial is returned through `d`. The same pattern as above can be used to reconstruct the expression. To extract the coefficient for a monomial term, indexing can be used. Of note, is an expression like `x^2*y*z` could *possibly* not equal the algebraically equal `x*y*z*x`, as they are only equal after some simplification, but the keys are in a canonical form, so this is not a concern:
```{julia}
d[x*y*z*x], d[z*y*x^2]
```
The residual term will capture any non-polynomial terms:
```{julia}
ex = sin(x) - x + x^3/6
d, r = polynomial_coeffs(ex, (x,))
r
```
To find the degree of a monomial expression, the `degree` function is available. Here it is applied to each monomial in `d`:
```{julia}
[degree(k) for (k,v) ∈ d]
```
The `degree` function will also identify the degree of more complicated terms:
```{julia}
degree(1 + x + x^2)
```
Mathematically the degree of the $0$ polynomial may be $-1$ or undefined, but here it is $0$:
```{julia}
degree(0), degree(1), degree(x), degree(x^a)
```
The coefficients are returned as *values* of a dictionary, and dictionaries are unsorted. To have a natural map between polynomials of a single symbol in the standard basis and a vector, we could use a pattern like this:
```{julia}
@variables x a0 as[1:10]
p = a0 + sum(as[i]*x^i for i ∈ eachindex(collect(as)))
d, r = polynomial_coeffs(p, (x,))
d
```
To sort the values we can use a pattern like the following:
```{julia}
vcat(r, [d[k] for k ∈ sort(collect(keys(d)), by=degree)])
```
---
Rational expressions can be decomposed into a numerator and denominator using the following idiom, which ensures the outer operation is division (a binary operation):
```{julia}
@variables x
ex = (1 + x + x^2) / (1 + x + x^2 + x^3)
function nd(ex)
ex1 = Symbolics.value(ex)
(operation(ex1) == /) || return (ex, one(ex))
Num.(arguments(ex1))
end
nd(ex)
```
With this, the study of asymptotic behaviour of a univariate rational expression would involve an investigation like the following:
```{julia}
m,n = degree.(nd(ex))
m > n ? "limit is infinite" : m < n ? "limit is 0" : "limit is a constant"
```
### Algebraically solving equations
The `~` operator creates a symbolic equation. For example
```{julia}
@variables x y
x^5 - x ~ 1
```
or
```{julia}
ex = [5x + 2y, 6x + 3y] .~ [1, 2]
```
The `Symbolics.solve_for` function can solve *linear* equations. For example,
```{julia}
Symbolics.solve_for(ex, [x, y])
```
### Limits
As of writing, there is no extra functionality provided by `Symbolics` for computing limits.
### Derivatives
`Symbolics` provides the `derivative` function to compute the derivative of a function with respect to a variable:
```{julia}
@variables a b c x
ex = a*x^2 + b*x + c
Symbolics.derivative(ex, x)
```
The computation can also be broken up into an expression indicating the derivative and then a function to apply the derivative rules:
```{julia}
D = Differential(x)
D(ex)
```
and then
```{julia}
expand_derivatives(D(ex))
```
The differentials can be multiplied to create operators for taking higher-order derivatives:
```{julia}
@variables x y
ex = (x - y^2)/(x^2 + y^2)
Dx, Dy = Differential(x), Differential(y)
Dxx, Dxy, Dyy = Dx*Dx, Dx*Dy, Dy*Dy
[Dxx(ex) Dxy(ex); Dxy(ex) Dyy(ex)] .|> expand_derivatives
```
In addition to `Symbolics.derivative` there are also the helper functions, such as `hessian` which performs the above
```{julia}
Symbolics.hessian(ex, [x,y])
```
The `gradient` function is also available
```{julia}
@variables x y z
ex = x^2 - 2x*y + z*y
Symbolics.gradient(ex, [x, y, z])
```
The `jacobian` takes an array of expressions:
```{julia}
@variables x y
exs = [ x^2 - y^2, 2x*y]
Symbolics.jacobian(exs, [x,y])
```
### Integration
The `SymbolicNumericIntegration` package provides a means to integrate *univariate* expressions through its `integrate` function.
Symbolic integration can be approached in different ways. SymPy implements part of the Risch algorithm in addition to other algorithms. Rules-based algorithms could also be implemented.
For example, here is a simple rule that could be used to integrate a single integral
```{julia}
is_var(x) = (xs = Symbolics.get_variables(x); length(xs) == 1 && xs[1] === x)
@syms x ∫(x)
r = @rule ∫(~x::is_var) => x^2/2
r(∫(x))
```
The `SymbolicNumericIntegration` package includes many more predicates for doing rules-based integration, but it primarily approaches the task in a different manner.
If $f(x)$ is to be integrated, a set of *candidate* answers is generated. The following is **proposed** as an answer: $\sum q_i \Theta_i(x)$. Differentiating the proposed answer leads to a *linear system of equations* that can be solved.
The example in the [paper](https://arxiv.org/pdf/2201.12468v2.pdf) describing the method is with $f(x) = x \sin(x)$ and the candidate thetas are ${x, \sin(x), \cos(x), x\sin(x), x\cos(x)}$ so that we propose:
$$
\int f(x) dx = q_1 x + q_2 \sin(x) + q_3 \cos(x) + q_4 x \sin(x) + q_4 x \cos(x)
$$
Differentiating both sides, yields a term $x\sin(x)$ on the left, and equating coefficients gives:
$$
q_1 = q_4 = 0,\quad q_5 = -1, \quad q_4 - q_3 = q_2 - q_5 = 0
$$
which can be solved with $q_5=-1$, $q_2=1$, and the other coefficients being $0$. That is $\int f(x) dx = 1 \sin(x) + (-1) x\cos(x)$.
The package provides an algorithm for the creation of candidates and the means to solve when possible. The `integrate` function is the main entry point. It returns three values: `solved`, `unsolved`, and `err`. The `unsolved` is the part of the integrand which can not be solved through this package. It is `0` for a given problem when `integrate` is successful in identifying an antiderivative, in which case `solved` is the answer. The value of `err` is a bound on the numerical error introduced by the algorithm.
To see, we have:
```{julia}
using SymbolicNumericIntegration
@variables x
integrate(x * sin(x))
```
The second term is `0`, as this has an identified antiderivative.
```{julia}
integrate(exp(x^2) + sin(x))
```
This returns `exp(x^2)` for the unsolved part, as this function has no simple antiderivative.
Powers of trig functions have antiderivatives, as can be deduced using integration by parts. When the fifth power is used, there is a numeric aspect to the algorithm that is seen:
```{julia}
u,v,w = integrate(sin(x)^5)
```
The derivative of `u` matches up to some numeric tolerance:
```{julia}
Symbolics.derivative(u, x) - sin(x)^5
```
The integration of rational functions (ratios of polynomials) can be done algorithmically, provided the underlying factorizations can be identified. The `SymbolicNumericIntegration` package has a function `factor_rational` that can identify factorizations.
```{julia}
import SymbolicNumericIntegration: factor_rational
@variables x
u = (1 + x + x^2)/ (x^2 -2x + 1)
v = factor_rational(u)
```
The summands in `v` are each integrable. We can see that `v` is a reexpression through
```{julia}
simplify(u - v)
```
The algorithm is numeric, not symbolic. This can be seen in these two factorizations:
```{julia}
u = 1 / expand((x^2-1)*(x-2)^2)
v = factor_rational(u)
```
or
```{julia}
u = 1 / expand((x^2+1)*(x-2)^2)
v = factor_rational(u)
```
As such, the integrals have numeric differences:
```{julia}
a,b,c = integrate(u)
```
We can see a bit of why through the following which needs a tolerance set to identify the rational numbers correctly:
```{julia}
cs = [first(arguments(term)) for term ∈ arguments(a)] # pick off coefficients
```
```{julia}
rationalize.(cs; tol=1e-8)
```

24
quarto/precalc/functions.qmd
View File

@@ -199,7 +199,7 @@ end
The conditions for the `if` statements are expressions that evaluate to either `true` or `false`, such as generated by the Boolean operators `<`, `<=`, `==`, `!-`, `>=`, and `>`.
If familiar with `if` conditions, they are natural to use. However, for simpler cases of "if-else" `Julia` provides the more convenient *ternary* operator: `cond ? if_true : if_false`. (The name comes from the fact that there are three arguments specified.) The ternary operator checks the condition and if true returns the first expression, whereas if the condition is false the second condition is returned. Both expressions are evaluated. (The [short-circuit](http://julia.readthedocs.org/en/latest/manual/control-flow/#short-circuit-evaluation) operators can be used to avoid both evaluations.)
If familiar with `if` conditions, they are natural to use. However, for simpler cases of "if-else" `Julia` provides the more convenient *ternary* operator: `cond ? if_true : if_false`. (The name comes from the fact that there are three arguments specified.) The ternary operator checks the condition and if true returns the first expression, whereas if the condition is false the second condition is returned. Both expressions are evaluated. (The [short-circuit](https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation) operators can be used to avoid both evaluations.)
For example, here is one way to define an absolute value function:
@@ -242,7 +242,7 @@ The `ternary` operator can be used to define an explicit domain. For example, a
```{julia}
hᵣ(t) = 0 <= t <= sqrt(10/16) ? 10.0 - 16t^2 : error("t is not in the domain")
h(t) = 0 <= t <= sqrt(10/16) ? 10.0 - 16t^2 : error("t is not in the domain")
```
#### Nesting ternary operators
@@ -252,7 +252,7 @@ The function `s(x)` isn't quite so easy to implement, as there isn't an "otherwi
```{julia}
s(x) = x < 0 ? 1 :
s(x) = x < 0 ? -1 :
x > 0 ? 1 : error("0 is not in the domain")
```
@@ -327,7 +327,7 @@ By using a multi-line function our work is much easier to look over for errors.
This next example, shows how using functions to collect a set of computations for simpler reuse can be very helpful.
An old method for finding a zero of an equation is the [secant method](https://en.wikipedia.org/wiki/Secant_method). We illustrate the method with the function $f(x) = x^2 - 2$. In an upcoming example we saw how to create a function to evaluate the secant line between $(a,f(a))$ and $(b, f(b))$ at any point. In this example, we define a function to compute the $x$ coordinate of where the secant line crosses the $x$ axis. This can be defined as follows:
An old method for finding a zero of an equation is the [secant method](https://en.wikipedia.org/wiki/Secant_method). We illustrate the method with the function $f(x) = x^2 - 2$. In an upcoming example we see how to create a function to evaluate the secant line between $(a,f(a))$ and $(b, f(b))$ at any point. In this example, we define a function to compute the $x$ coordinate of where the secant line crosses the $x$ axis. This can be defined as follows:
```{julia}
@@ -470,10 +470,10 @@ mxb(0)
So the `b` is found from the currently stored value. This fact can be exploited. we can write template-like functions, such as `f(x)=m*x+b` and reuse them just by updating the parameters separately.
How `Julia` resolves what a variable refers to is described in detail in the manual page [Scope of Variables](http://julia.readthedocs.org/en/latest/manual/variables-and-scoping/). In this case, the function definition finds variables in the context of where the function was defined, the main workspace. As seen, this context can be modified after the function definition and prior to the function call. It is only when `b` is needed, that the context is consulted, so the most recent binding is retrieved. Contexts (more formally known as environments) allow the user to repurpose variable names without there being name collision. For example, we typically use `x` as a function argument, and different contexts allow this `x` to refer to different values.
How `Julia` resolves what a variable refers to is described in detail in the manual page [Scope of Variables](https://docs.julialang.org/en/v1/manual/variables-and-scoping/). In this case, the function definition finds variables in the context of where the function was defined, the main workspace. As seen, this context can be modified after the function definition and prior to the function call. It is only when `b` is needed, that the context is consulted, so the most recent binding is retrieved. Contexts (more formally known as environments) allow the user to repurpose variable names without there being name collision. For example, we typically use `x` as a function argument, and different contexts allow this `x` to refer to different values.
Mostly this works as expected, but at times it can be complicated to reason about. In our example, definitions of the parameters can be forgotten, or the same variable name may have been used for some other purpose. The potential issue is with the parameters, the value for `x` is straightforward, as it is passed into the function. However, we can also pass the parameters, such as $m$ and $b$, as arguments. For parameters, we suggest using [keyword](http://julia.readthedocs.org/en/latest/manual/functions/#keyword-arguments) arguments. These allow the specification of parameters, but also give a default value. This can make usage explicit, yet still convenient. For example, here is an alternate way of defining a line with parameters `m` and `b`:
Mostly this works as expected, but at times it can be complicated to reason about. In our example, definitions of the parameters can be forgotten, or the same variable name may have been used for some other purpose. The potential issue is with the parameters, the value for `x` is straightforward, as it is passed into the function. However, we can also pass the parameters, such as $m$ and $b$, as arguments. For parameters, we suggest using [keyword](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments) arguments. These allow the specification of parameters, but also give a default value. This can make usage explicit, yet still convenient. For example, here is an alternate way of defining a line with parameters `m` and `b`:
```{julia}
@@ -527,7 +527,7 @@ For example, here we use a *named tuple* to pass parameters to `f`:
```{julia}
#| hold: true
function trajectory(x ,p)
g,v0, theta, k = p.g, p.v0, p.theta, p.k # unpack parameters
g, v0, theta, k = p.g, p.v0, p.theta, p.k # unpack parameters
a = v0 * cos(theta)
(g/(k*a) + tan(theta))* x + (g/k^2) * log(1 - k/a*x)
@@ -663,7 +663,7 @@ With `Julia` we can represent such operations. The simplest thing would be to do
```{julia}
#| hold: true
f(x) = x^2 - 2x
g(x) = f(x -3)
g(x) = f(x - 3)
```
Then $g$ has the graph of $f$ shifted by 3 units to the right. Now `f` above refers to something in the main workspace, in this example a specific function. Better would be to allow `f` to be an argument of a function, like this:
@@ -731,7 +731,7 @@ To model this in `Julia`, we would want to turn the inputs `f`,`a`, `b` into a f
```{julia}
function secant(f, a, b)
m = (f(b) - f(a)) / (b-a)
m = (f(b) - f(a)) / (b - a)
x -> f(a) + m * (x - a)
end
```
@@ -767,7 +767,7 @@ specific_line(m,b) = x -> mxplusb(x; m=m, b=b)
The returned object will have its parameters (`m` and `b`) fixed when used.
In `Julia`, the functions `Base.Fix1` and `Base.Fix2` are provided to take functions of two variables and create callable objects of just one variable, with the other argument fixed. This partial function application is provided by a some of the logical comparison operators. which can be useful with filtering, say.
In `Julia`, the functions `Base.Fix1` and `Base.Fix2` are provided to take functions of two variables and create callable objects of just one variable, with the other argument fixed. This partial function application is provided by a some of the logical comparison operators, which can be useful with filtering, say.
For example, `<(2)` is a funny looking way of expressing the function `x -> x < 2`. (Think of `x < y` as `<(x,y)` and then "fix" the value of `y` to be `2`.) This is useful with filtering by a predicate function, for example:
@@ -1010,7 +1010,7 @@ A(w, h) = w * h
when called as `A(10, 5)` will use 10 for `w` and `5` for `h`, as the order of `w` and `h` matches that of `10` and `5` in the call.
This is clear enough, but in fact positional arguments can have default values (then called [optional](http://julia.readthedocs.org/en/latest/manual/functions/#optional-arguments)) arguments). For example,
This is clear enough, but in fact positional arguments can have default values (then called [optional](https://docs.julialang.org/en/v1/manual/functions/#Optional-Arguments)) arguments). For example,
```{julia}
@@ -1132,7 +1132,7 @@ radioq(choices, answ, keep_order=true)
###### Question
The [pipe](http://julia.readthedocs.org/en/latest/stdlib/base/#Base.|>) notation `ex |> f` takes the output of `ex` and uses it as the input to the function `f`. That is composition. What is the value of this expression `1 |> sin |> cos`?
The [pipe](https://docs.julialang.org/en/v1/manual/functions/#Function-composition-and-piping) notation `ex |> f` takes the output of `ex` and uses it as the input to the function `f`. That is composition. What is the value of this expression `1 |> sin |> cos`?
```{julia}

2
quarto/precalc/logical_expressions.qmd
View File

@@ -47,7 +47,7 @@ The set of numeric comparisons is nearly the same as the mathematical counterpar
:::{.callout-warning}
## Warning
The use of `==` is necessary, as `=` is used for assignment and mutation.")
The use of `==` is necessary, as `=` is used for assignment and mutation.
:::

10
quarto/precalc/numbers_types.qmd
View File

@@ -58,8 +58,8 @@ Though there are explicit constructors for these types, these notes avoid them u
:::{.callout-note}
## Warngin
Heads up, the difference between `1` and `1.0` is subtle. Even more so, as `1.` will parse as `1.0`. This means some expressions, such as `2.*3`, are ambigous, as the `.` might be part of the `2` (as in `2. * 3`) or the operation `*` (as in `2 .* 3`).
## Warning
Heads up, the difference between `1` and `1.0` is subtle. Even more so, as `1.` will parse as `1.0`. This means some expressions, such as `2.*3`, are ambiguous, as the `.` might be part of the `2` (as in `2. * 3`) or the operation `*` (as in `2 .* 3`).
:::
@@ -69,7 +69,7 @@ Similarly, each type is printed slightly differently.
The key distinction is between integers and floating points. While floating point values include integers, and so can be used exclusively on the calculator, the difference is that an integer is guaranteed to be an exact value, whereas a floating point value, while often an exact representation of a number is also often just an *approximate* value. This can be an advantage – floating point values can model a much wider range of numbers.
Now in nearly all cases the differences are not noticable. Take for instance this simple calculation involving mixed types.
Now in nearly all cases the differences are not noticeable. Take for instance this simple calculation involving mixed types.
```{julia}
@@ -200,14 +200,14 @@ or may be represented in scientific notation:
```{julia}
6.23 * 10.0^23
6.022 * 10.0^23
```
The special coding `aeb` (or if the exponent is negative `ae-b`) is used to represent the number $a \cdot 10^b$ ($1 \leq a < 10$). This notation can be used directly to specify a floating point value:
```{julia}
avagadro = 6.23e23
avogadro = 6.022e23
```
Here `e` is decidedly *not* the Euler number, rather syntax to separate the exponent from the mantissa.

10
quarto/precalc/vectors.qmd
View File

@@ -410,13 +410,11 @@ The last value of a vector is usually denoted by $v_n$. In `Julia`, the `length`
There is [much more](http://julia.readthedocs.org/en/latest/manual/arrays/#indexing) to indexing than just indexing by a single integer value. For example, the following can be used for indexing:
* a scalar integer (as seen)
:::
* a range
* a vector of integers
* a boolean vector
:::
Some add-on packages extend this further.
@@ -541,7 +539,7 @@ sum(v), length(v)
Other desired operations with vectors act differently. Rather than reduce a collection of values using some formula, the goal is to apply some formula to *each* of the values, returning a modified vector. A simple example might be to square each element, or subtract the average value from each element. An example comes from statistics. When computing a variance, we start with data $x_1, x_2, \dots, x_n$ and along the way form the values $(x_1-\bar{x})^2, (x_2-\bar{x})^2, \dots, (x_n-\bar{x})^2$.
Such things can be done in *many* differents ways. Here we describe two, but will primarily utilize the first.
Such things can be done in *many* different ways. Here we describe two, but will primarily utilize the first.
### Broadcasting a function call
@@ -678,7 +676,7 @@ The first task is to create the data. We will soon see more convenient ways to g
```{julia}
a,b, n = -1, 1, 7
a, b, n = -1, 1, 7
d = (b-a) // (n-1)
𝐱s = [a, a+d, a+2d, a+3d, a+4d, a+5d, a+6d] # 7 points
```
@@ -712,7 +710,7 @@ The style generally employed here is to use plural variable names for a collecti
## Other container types
Vectors in `Julia` are a container, one of many different types. Another useful type for programming purposes are *tuples*. If a vector is formed by placing comma-separated values within a `[]` pair (e.g., `[1,2,3]`), a tuple is formed by placing comma-separated values withing a `()` pair. A tuple of length $1$ uses a convention of a trailing comma to distinguish it from a parethesized expression (e.g. `(1,)` is a tuple, `(1)` is just the value `1`).
Vectors in `Julia` are a container, one of many different types. Another useful type for programming purposes are *tuples*. If a vector is formed by placing comma-separated values within a `[]` pair (e.g., `[1,2,3]`), a tuple is formed by placing comma-separated values withing a `()` pair. A tuple of length $1$ uses a convention of a trailing comma to distinguish it from a parenthesized expression (e.g. `(1,)` is a tuple, `(1)` is just the value `1`).
Tuples are used in programming, as they don't typically require allocated memory to be used so they can be faster. Internal usages are for function arguments and function return types. Unlike vectors, tuples can be heterogeneous collections. (When commas are used to combine more than one output into a cell, a tuple is being used.) (Also, a big technical distinction is that tuples are also different from vectors and other containers in that tuple types are *covariant* in their parameters, not *invariant*.)

1
.github/workflows/publish.yml → quarto/publish.yml
View File

@@ -1,3 +1,4 @@
# This is *if* we have CI publish; currently we use `quarto render gh-pages`
on:
push:
branches: main