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

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.

Comparativa de la sintaxis de Julia con otros lenguajes

La sintaxis de Julia es similar a la de lenguajes como Python o MATLAB, y en algunos casos es incluso más concisa, por lo que la transición de estos lenguajes a Julia es fácil.

https://cheatsheets.quantecon.org

Algunas razones más para pensárselo

Evita el problema de la doble reescritura de código (prototipo + versión operativa).
Permite la programación paralela y distribuida.
Puede ejecutar código realizado en otros lenguajes como C, Fortran, Python, R, Matlab, etc.

using PyCall, RCall 
@pyimport scipy.interpolate as interpolate

xs = 1:10
ys = rand(10)
f = interpolate.interp1d(xs, ys, kind="cubic")
zs = range(1, 10, 1000)

R"""
plot($xs, $ys)
lines($zs, $(f(zs)), col = "red")
"""

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)

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

😄 = "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{PyCall.PyObject}

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\)

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))

Set{Sym} with 2 elements:
  0
  1

solveset(x^2 - x)

Set{Sym} with 2 elements:
  0
  1

solveset(sin(x) - 1)

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

solveset(x^2 - y, x)

Set{Sym} with 2 elements:
  -sqrt(y)
  sqrt(y)

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

Set{Sym} with 1 element:
  (2, 1)

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
SymPy.limit(sin(x) / x, x=>0)

\(1\)

SymPy.limit(1 / x, x=>0, dir = "-")

\(-\infty\)

SymPy.limit(1 / x, x=>0, dir = "+")

\(\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} e^{x y} + 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))

\(\left[\begin{smallmatrix}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{smallmatrix}\right]\)

v = [diff.(f, [x, y]) for f in [x^2*y, x*y^2]]
hcat(v...)

\(\left[\begin{smallmatrix}2 x y & y^{2}\\x^{2} & 2 x y\end{smallmatrix}\right]\)

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} & e^{x y} + x y e^{x y} \\ e^{x y} + 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, 1.1102230246251565e-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}\)