Cálculo numérico y simbólico con Julia

Seminario de Innovación Docente
Universidade da Coruña

Alfredo Sánchez Alberca  

21 de julio de 2023

¿Qué es Julia?

Julia un lenguaje de programación moderno creado en 2012 el MIT por el equipo del profesor Edelman, orientado a cálculo científico y el análisis de datos.

De los creadores de Julia:

We want a language that is:

  • Open source.
  • With the speed of C.
  • Obvious, familiar mathematical notation like Matlab.
  • As usable for general programming as Python.
  • As easy for statistics as R.
  • As natural for string processing as Perl.
  • As powerful for linear algebra as Matlab.
  • As good at gluing programs together as the shell.
  • Dirt simple to learn, yet keeps the most serious hackers happy.

Experiencia docente

Prácticas en el Grado de Ingeniería Matemática

Curso 2022/2023

  • Álgebra Lineal
  • Análisis I y II
  • Matemática Discreta

Curso 2023/2024

  • Análisis III
  • Ecuaciones Diferenciales
  • Ecuaciones en Derivadas Parciales
  • Geometría Diferencial
  • Métodos numéricos I y II

Curso 2024/2025

  • Aprendizaje automático
  • Redes neuronales
  • Series temporales
  • Procesamiento del lenguaje natural

¿Por qué Julia?

Sencillez

Lenguaje de alto nivel con una sintaxis fácil de aprender (similar a Python, R o Matlab).

println("¡Hola Julia!")
¡Hola Julia!
x = [1, 2, 3, 4, 5]
n = length(x)  
media = sum(x) / n
varianza = sum(x.^2) / n - media^2
(media, varianza)
(3.0, 2.0)

Rapidez

Julia es un lenguaje muy veloz (equiparable a C o Fortran en muchas tareas).

Comparativa de Julia con otros lenguajes

Cálculo de pi por el método de Leibniz

Julia

function pi(n)
  pi = 0.0;
  op = 1
  for i in 0:n
      pi += op*4/(2i+1)
      op *= -1
  end
  return pi  
end

pi(10^9)

\(\approx\) 1 seg

C

double pi(double n) {
  double pi = 0;
  int i = 0;
  int op = 1;
  for (i; i < n; i++) {
    pi += op*4/(2*i+1);
    op *= -1;
  }
  return pi;
}

int main() {
  pi(1000000000);
}

\(\approx\) 3 seg

Python

def pi(n): 
    pi = 0;
    op = 1
    for i in range(n):
        pi += op*4/(2*i+1)
        op *= -1
    return pi 

pi(10^9)

\(\approx\) 78 seg

Despacho múltiple

La razón por la que Julia es tan rápido es por que usa despacho múltiple para precompilar distintas versiones de una misma función para cada posible tipo de dato de sus argumentos.

Tipos de datos numéricos

typeof(1)
Int64
typeof(Int16(1))
Int16
typeof(1.0)
Float64
typeof(1.0f0)
Float32
typeof(1//2)
Rational{Int64}
typeof(pi)
Irrational{:π}
typeof(1+1im)
Complex{Int64}

Enteros

Tipo Signo Bits Min Max
Int8 ✓ 8 \(-2^7\) \(2^7-1\)
UInt8 8 \(0\) \(2^8-1\)
Int16 ✓ 16 \(-2^{15}\) \(2^{15}-1\)
UInt16 16 \(0\) \(2^{16}-1\)
Int32 ✓ 32 \(-2^{31}\) \(2^{31}-1\)
UInt32 32 \(0\) \(2^{32}-1\)
Int64 ✓ 64 \(-2^{63}\) \(2^{63}-1\)
UInt64 64 \(0\) \(2^{64}-1\)
Int128 ✓ 128 \(-2^{127}\) \(2^{127}-1\)
UInt128 128 \(0\) \(2^{128}-1\)
Bool N/A 8 false (0) true (1)

Reales en coma flotante

Tipo Bits Min Max
Float16 16 \(-6.55e4\) \(6.55e4\)
Float32 32 \(-3.4028235e38\) \(3.4028235e38\)
Float64 64 \(-1.7976931348623157e308\) \(1.7976931348623157e308\)
Tipo Min positive Epsilon
Float16 \(6.104e-5\) \(0.000977\)
Float32 \(1.1754944f-38\) \(1.1920929f-7\)
Float64 \(2.2250738585072014e-308\) \(2.220446049250313e-16\)

Precisión arbitraria

  • BigInt permite representar enteros de cualquier tamaño.
  • BigFloat permite representar reales de cualquier tamaño y precisión.

El límite lo marca la memoria del ordenador.

BigFloat(pi)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198
BigFloat(pi, precision = 1000)
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412736

Definición de nuevos tipos

struct EvenInt <: Integer
    num :: Signed

    # inner constructor
    function EvenInt(num)
        @assert iseven(num)
        new(num)
    end
end

EvenInt(4)
EvenInt(4)

Aritmética racional

  • Simplifiación

    4 // 6
    2//3
    4 // 6 == 2 // 3
    true

  • Suma y resta

    1 // 6 + 2 // 4
    2//3
    2 // 3 - 1 // 6
    1//2

  • Producto

    2 // 3 * 3 // 4
    1//2

  • Cociente

    3 // 4 / 2 // 3
    9//8

Aritmética compleja

La constante im está asociada al número complejo \(i=\sqrt{-1}\).

  • Suma y resta

    (1 + 2im) + (1 - im)
    2 + 1im
    (1 + 2im) - (1 + im)
    0 + 1im

  • Producto

    (1 + 2im) * (2 - 3im)
    8 + 1im

  • Cociente

    (1 + 2im) / (1 - 2im)
    -0.6 + 0.8im

Promoción de tipos

Siempre se promociona hacia el tipo más general.

typeof(Int16(1) + 2)
Int64
typeof(1.0 + 2)
Float64
typeof(1 / 2)
Float64
typeof(1 // 2 + 0.5)
Float64
typeof(sqrt(1 // 2))
Float64
typeof((1 // 2) ^ 2.0)
Float64
typeof(1.0 + 2)
Float64
typeof(2(1+0im))
Complex{Int64}
typeof(im^2)
Complex{Int64}

Tipos de datos compuestos

  • Vectores

    typeof([1, 2, 3])
    Vector{Int64} (alias for Array{Int64, 1})

  • Matrices

    typeof([1 2 3; 4 5 6])
    Matrix{Int64} (alias for Array{Int64, 2})

  • Tuplas

    typeof((1, "enero", 2020))
    Tuple{Int64, String, Int64}

  • Diccionarios

    typeof(Dict("x" => 1, "y" => 2))
    Dict{String, Int64}

  • Conjuntos

    typeof(Set([2, 4, 6]))
    Set{Int64}

Los tipos compuestos incluyen parámetros que determinan el tipo de sus elementos.

Álgebra de conjuntos

A = Set([2, 4])
B = Set([1, 2, 3])
3 ∈ A 
false
A ⊆ B
false
A ∪ B
Set{Int64} with 4 elements:
  4
  2
  3
  1
A ∩ B
Set{Int64} with 1 element:
  2

Es fácil definir nuevos operadores.

const \ˢ = setdiff
A \ˢ B
Set{Int64} with 1 element:
  4

Aritmética infinita e indeterminaciones

1 / Inf
0.0
1 / 0
Inf
0 / 0
NaN
Inf + 1
Inf
1 - Inf
-Inf
Inf + Inf
Inf
Inf - Inf
NaN
Inf / Inf
NaN
0 / 0
NaN
0 * Inf
NaN

Definición de funciones

Julia admite múltiples formas de definir funciones.

  • Tradicional. Similar a otros lenguajes de programación.

    function area_triángulo(base, altura)
  base * altura / 2
end

  • En línea. Muy cómoda para definir funciones matemáticas.

    area_triángulo(base, altura) = base * altura / 2

  • Funciones anónimas. Se utilizan sobre todo en programación funcional, como argumentos de otras funciones.

    area_triángulo = (base, altura) -> base * altura / 2

Programación funcional

Las funciones son objetos que pueden asignarse o pasarse como argumentos de otras funciones.

square(x) = x * x
cuadrado = square
cuadrado(2)
4
componer(f, g) = f ∘ g
componer(log, sin)(pi/2)
0.0

Vectorización de funciones y operadores

El operador punto . permite aplicar una función u operador a los elementos de una colección (arrays, vectores, tuplas, etc.)

v = [1, 2, 3]
cuadrado.(v)
3-element Vector{Int64}:
 1
 4
 9
w = [0, 1, 2]
v .^ w
3-element Vector{Int64}:
 1
 2
 9

Expresividad matemática

Acepta Unicode

Julia permite caracteres Unicode lo que facilita la expresión de fórmulas matemáticas.

Para ello se utilizan códigos especiales (en muchos casos son los mismos que en \(\LaTeX\)), pulsando después la tecla de tabulación.

# Introducir \alpha + TAB
α = 1
1
😄 = "julia"
😄 ^ 2
"juliajulia"

Composición de funciones

El operador ∘ permite componer funciones.

f(x) = sin(x)
g(x) = log(x)
(g ∘ f)(pi/2)
0.0

También puede realizarse con el operador de tubería |> (similar a %>% en R).

pi / 2 |> sin |> log
0.0

Integración con \(\LaTeX\)

El paquete Latexify permite obtener el código \(\LaTeX\) para representar cualquier expresión matemática.

using Latexify
expr = :((x-y)/(x+y)^2)
latexify(expr)

\(\frac{x - y}{\left( x + y \right)^{2}}\)

Integración con Quarto

Quarto es un sistema de publicación de textos científicos y técnicos de alta calidad basado en Markdown y \(\LaTeX\) (similar a RMarkdown).

Julia es uno de los lenguajes de programación que admite.

Orientación al cálculo numérico

Lenguaje de propósito general, pero especialmente diseñado para el cálculo científico y el análisis de datos.

Aplicaciones Matemáticas

Representación gráfica de funciones

El paquete Plots permite la representación gráfica de funciones con diversos tipos de diagramas.

using Plots
plot(sin, 0, 2pi, label = "sen", size = (480,400))
plot!(cos, label = "cos")
xs = ys = range(-3,3, length=100)
f(x,y) = 2sin(x^2+y^2) / (x^2+y^2)
surface(xs, ys, f, size = (480, 400))
plotlyjs()
xs = ys = range(-3,3, length=100)
f(x,y) = 2sin(x^2+y^2) / (x^2+y^2)
surface(xs, ys, f)

The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.

Cálculo simbólico

Existen dos paquetes para cálculo simbólico:

  • SymPy: Basado en el paquete sympy de Phthon. Muy maduro, pero lento.
  • Symbolics. Paquete puro en Julia muy rápido, pero no tan completo.

Declaración de variables (símbolos)

SymPy

using SymPy
@syms x::real y::real
typeof(x)
Sym

Symbolics

using Symbolics
@variables x y;
typeof(x)
Num

Simplificación

SymPy

using SymPy
@syms x::real
SymPy.simplify(sin(x)^2 + cos(x)^2)

\(1\)

using SymPy
@syms x::real
SymPy.expand((x + 1)^2)

\(x^{2} + 2 x + 1\)

SymPy.factor(x^3 - x^2 + x - 1)

\(\left(x - 1\right) \left(x^{2} + 1\right)\)

Symbolics

using Symbolics
@variables x;
Symbolics.simplify(sin(x)^2 + cos(x)^2)

\[ \begin{equation} 1 \end{equation} \]

Sustitución de variables

SymPy

using SymPy
@syms x y
expr = cos(x*y)
expr.subs([(x, 2), (y, π)])

\(1.0\)

expr(x=>2, y=>π)

\(1\)

Symbolics

using Symbolics
@variables x y
expr = cos(x*y)
substitute(expr, Dict([x=>2, y=>π]))

\[ \begin{equation} 1 \end{equation} \]

Resolución de ecuaciones

SymPy

using SymPy
@syms x y
solveset(Eq(x^2, x))

\(\left\{0, 1\right\}\)

solveset(x^2 - x)

\(\left\{0, 1\right\}\)

solveset(sin(x) - 1)

\(\left\{2 n \pi + \frac{\pi}{2}\; \middle|\; n \in \mathbb{Z}\right\}\)

solveset(x^2 - y, x)

\(\left\{- \sqrt{y}, \sqrt{y}\right\}\)

linsolve((x+y-3, x-y-1), (x,y))

\(\left\{\left( 2, \ 1\right)\right\}\)

Symbolics
Para sistemas lineales

using Symbolics
@variables x y
Symbolics.solve_for([x + y ~ 3, x - y ~ 1], [x, y])
2-element Vector{Float64}:
 2.0
 1.0

Límites

SymPy

using SymPy
@syms x
limit(sin(x) / x, x, 0)

\(1\)

limit(sin(x) / x, x=>0)

\(1\)

limit(1 / x, x, 0, "-")

\(-\infty\)

limit(1 / x, x, 0, "+")

\(\infty\)

Symbolics
No implementado aún.

Derivadas

SymPy

using SymPy
@syms x y
diff(sin(x))

\(\cos{\left(x \right)}\)

diff(sin(x), x, 2)

\(- \sin{\left(x \right)}\)

SymPy.diff(ℯ^(x*y), x)

\(y e^{x y}\)

diff(ℯ^(x*y), x, y)

\(\left(x y + 1\right) e^{x y}\)

Symbolics

using Symbolics
@variables x y
Symbolics.derivative(sin(x), x)

\[ \begin{equation} \cos\left( x \right) \end{equation} \]

Dx = Symbolics.Differential(x)
expand_derivatives(Dx(sin(x)))

\[ \begin{equation} \cos\left( x \right) \end{equation} \]

Dy = Symbolics.Differential(y)
expand_derivatives(Dy(Dx(ℯ^(x*y))))

\[ \begin{equation} x y e^{x y} + e^{x y} \end{equation} \]

Es posible usar también la notación prima f' sobrecargando el operador adjoin.

Gradiente, Hessiana y Jacobiana

SymPy

using SymPy
@syms x y
diff.(ℯ^(x*y), (x, y))
(y*exp(x*y), x*exp(x*y))
hessian(ℯ^(x*y), (x,y))
2×2 Matrix{Sym}:
            y^2*exp(x*y)  x*y*exp(x*y) + exp(x*y)
 x*y*exp(x*y) + exp(x*y)             x^2*exp(x*y)
v = [diff.(f, [x, y]) for f in [x^2*y, x*y^2]]
hcat(v...)
2×2 Matrix{Sym}:
 2⋅x⋅y    y^2
   x^2  2⋅x⋅y

Symbolics

using Symbolics
@variables x y
Symbolics.gradient(ℯ^(x*y), [x, y])

\[ \begin{equation} \left[ \begin{array}{c} y e^{x y} \\ x e^{x y} \\ \end{array} \right] \end{equation} \]

Symbolics.hessian(ℯ^(x*y), [x, y])

\[ \begin{equation} \left[ \begin{array}{cc} y^{2} e^{x y} & x y e^{x y} + e^{x y} \\ x y e^{x y} + e^{x y} & x^{2} e^{x y} \\ \end{array} \right] \end{equation} \]

Symbolics.jacobian([ℯ^(x*y), x*y], [x, y])

\[ \begin{equation} \left[ \begin{array}{cc} y e^{x y} & x e^{x y} \\ y & x \\ \end{array} \right] \end{equation} \]

Integrales

SymPy

using SymPy
@syms x y
integrate(cos(x))

\(\sin{\left(x \right)}\)

integrate(cos(x), (x, 0, π/2))

\(1.0\)

integrate(x^2 - y, (x, -1, 1), (y, -1, 1))

\(\frac{4}{3}\)

Symbolics

using Symbolics
using SymbolicNumericIntegration
@variables x y
SymbolicNumericIntegration.integrate(cos(x))
(sin(x), 0, 0)
SymbolicNumericIntegration.integrate(x^2 / (16 + x^2))
(x - 4atan((1//4)*x), 0, 2.2215299868541707e-16)

Ecuaciones diferenciales

using SymPy
@syms x k
y = SymFunction("y")
dsolve(diff(y(x),x) ~  k*y(x))

\(y{\left(x \right)} = C_{1} e^{k x}\)

dsolve(diff(y(x),x) ~  k*y(x), ics = Dict(y(0) => 1))

\(y{\left(x \right)} = e^{k x}\)

Cálculo numérico

Polinomios

using Polynomials
p = Polynomial([4,0,2,1])
p(2)
20
Polynomials.roots(p)
3-element Vector{ComplexF64}:
 -2.5943130163548487 + 0.0im
 0.29715650817742467 - 1.205625150602913im
 0.29715650817742467 + 1.205625150602913im
fit([1,2,3], [1, 0, 1])
4.0 - 4.0∙x + 1.0∙x2

Interpolación

using Interpolations
xs = 1:10
itp = Interpolations.interpolate(log.(xs), BSpline(Linear()))
10-element interpolate(::Vector{Float64}, BSpline(Linear())) with element type Float64:
 0.0
 0.6931471805599453
 1.0986122886681098
 1.3862943611198906
 1.6094379124341003
 1.791759469228055
 1.9459101490553132
 2.0794415416798357
 2.1972245773362196
 2.302585092994046
xs = 1:0.5:10
itp2 = Interpolations.interpolate(log.(xs), BSpline(Quadratic(Line(OnGrid()))))
19-element interpolate(OffsetArray(::Vector{Float64}, 0:20), BSpline(Quadratic(Line(OnGrid())))) with element type Float64:
 0.0
 0.4054651081081644
 0.6931471805599454
 0.9162907318741551
 1.0986122886681098
 1.252762968495368
 1.3862943611198908
 1.504077396776274
 1.6094379124341005
 1.7047480922384253
 1.791759469228055
 1.8718021769015913
 1.9459101490553135
 2.0149030205422642
 2.079441541679836
 2.1400661634962708
 2.1972245773362196
 2.2512917986064953
 2.302585092994046

Álgebra Lineal

A = [1 2 3; 0 1 0; 1 0 1]
B = [10, 2, 4]
A \ B
3-element Vector{Float64}:
 3.0
 2.0
 1.0
using LinearAlgebra
eigvals(A)
3-element Vector{Float64}:
 -0.7320508075688774
  1.0
  2.732050807568877
factorize(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0   0.0  0.0
 1.0   1.0  0.0
 0.0  -0.5  1.0
U factor:
3×3 Matrix{Float64}:
 1.0   2.0   3.0
 0.0  -2.0  -2.0
 0.0   0.0  -1.0

Resolución de ecuaciones

using Roots
f(x) = sin(x)
find_zeros(f, -5, 5)
3-element Vector{Float64}:
 -3.141592653589793
  0.0
  3.141592653589793
using LinearSolve
A = [1.0 2 3; 0 1 0; 1 0 1]
B = [10.0, 2, 4]
solve(LinearProblem(A, B))
u: 3-element Vector{Float64}:
 3.0
 2.0
 1.0
using NLsolve
f(x, y) = [(x+3)*(y^3-7)+18, sin(y*exp(x)-1)]
f(x) = [(x[1]+3)*(x[2]^3-7)+18, sin(x[2]*exp(x[1])-1)]
nlsolve(f, [ 0.1, 1.2])
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.1, 1.2]
 * Zero: [-7.775548712324193e-17, 0.9999999999999999]
 * Inf-norm of residuals: 0.000000
 * Iterations: 4
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 5
 * Jacobian Calls (df/dx): 5

Derivación numérica

using ForwardDiff
f(x) = exp(-x^2)
ForwardDiff.derivative(f, 1)
-0.7357588823428847

Integración numérica

using QuadGK
f(x) = x^x
quadgk(f, 0, 1)
(0.7834305106229383, 4.864047812094818e-9)
quadgk(f, big(0), big(1), rtol = 1e-30)
(0.7834305107121344070592643865269751743063670117483685034318247197232693913877173, 7.736858106873745132469365322882480666696649544982787019412483807349385217202092e-31)
quadgk(x -> ℯ^(-x^2), -Inf, Inf)
(1.7724538509055137, 6.4296367091810765e-9)

Ecuaciones diferenciales

using OrdinaryDiffEq
S0, I0, R0 = 7900000, 10, 0
N = S0 + I0 + R0
u0 = [S0, I0, R0]/N
k = 1/3
f(u,p,t) = -k * u  # solving  u′(t) = - k u(t)
time_span = (0.0, 20.0)
prob = ODEProblem(f, I0/N, time_span)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plot(sol)