L'idea è che ad ogni articoletto di "divulgazione matematica" dovrebbe essere associato un video che lo racconta/analizza/risolve in dettaglio. Se non c'è è perché devo ancora registrarlo :)
Presto inoltre ristrutturerò la pagina, il suo layout, in modo che ogni problema abbia la sua pagina dedicata, per poi scriverci meglio tutto all'interno.
Finché si lavora con quadrati o cubi è tutto molto easy, possiamo anche disegnarli per studiarli e per capire, ad esempio, quanti vertici, lati e facce abbiano. Salendo di dimensione occorre invece arrivarci per vie traverse, e qui vediamo come ci si può arrivare sfruttando anche un po' di calcolo combinatorio.
Queste sarebbero le conclusioni a cui si arriva: detto \(N\) il numero di dimensioni in cui lavoriamo (quindi considerando un \(N\)-cubo), abbiamo che
\[ \begin{align*} \#\text{(vertici)} &= 2^{N}\\ \#\text{(lati)} &= \dfrac{\#\text{(vertici)}\cdot N}{2}\\ \#\text{(facce piane)} &= \dfrac{\#\text{(vertici)}\cdot \binom{N}{2}}{2^2}\\ \#\text{(facce di dimensione $k$)} &= \dfrac{\#\text{(vertici)}\cdot \binom{N}{k}}{2^k} \quad k\in\{3,\ldots,N-1\} \end{align*} \]dove \(\binom{N}{k}\) è il coefficiente binomiale di \(N\) su \(k\), quindi \(\frac{N!}{k!(N-k)!}\). Le facce piane sono quelle classiche, mentre quelle di dimensioni maggiori hanno senso in dimensioni maggiori. Per esempio, uno può chiedersi da quante facce 5-cubiche è formato un 7-cubo, o domande ancora peggiori.
Questo problema è noto come "paradosso dei compleanni" perché per la prima probabilità si arriva al 50% già con solo 23 persone, che in effetti sembrano molte poche. Nel video è quindi spiegato come ricavare quei valori e anche come interpretarli, per capire perché quel risultato un po' controintuitivo abbia in realtà senso.
Sempre nel caso della prima probabilità si raggiuge davvero l'1 (ovvero la certezza che si realizzi l'evento) ovviamente solo con 366 persone (per il principio della piccionaia, volendo essere precisi, in modo che vengono forzate le date di due persone ad essere uguali). Gli "uno" prima di lui presenti nel plot che ora segue sono quindi solo dovuti ad approssimazioni del grafico (che non riporta tutti i decimali) che diventano poi approssimazioni numeriche (in cui cioè i calcoli saturano la precisione disponibile per la rappresentazione dei numeri decimali in un computer) da 295 persone in poi.
Infatti i valori precisi calcolati dal codice sarebbero:
⋮
n=80: 0.9999143319493134946903219036781801337436248827650357978945842589451803561427634
n=81: 0.9999331085083680711965527193103598304573509358576306915067301747928120589059932
n=82: 0.9999479529215795403282766363949101146572264815988140175011270401127633554227453
⋮
n=258: 0.9999999999999999999999999999999999999999999999999999998264153107735258302391866
n=259: 0.9999999999999999999999999999999999999999999999999999999491135294596363392755936
n=260: 0.999999999999999999999999999999999999999999999999999999985222011295127265652638
⋮
n=293: 0.9999999999999999999999999999999999999999999999999999999999999999999999999999309
n=294: 0.9999999999999999999999999999999999999999999999999999999999999999999999999999827
n=295: 1.0
⋮using Plots
p1(n) = 1 - factorial(big(365)) / (factorial(big(365-n)) * big(365)^n)
p1_values = p1.(1:365)
plot(p1_values,title="p1 values")
p2(n) = 1 - (big(364/365)^n)
p2_values = p2.(1:2000)
plot(p2_values,title="p2 values")Vediamo come calcolare la probabilità di vincere un tiebreak a tennis, quindi una super applicazione della matematica (statistica e probabilità) ad un problema reale.
using Plots
using Distributions
# 1 2 3 4 5 6 7 8 9 10 11 12
# A B B A A B B A A B B A
function ptiebreak(a::Real,b::Real)
p70 = (pdf(Binomial(3,a),3) * pdf(Binomial(3,b),3))*b
# batte B al 7mo turno
p71 = sum([pdf(Binomial(3,a),i) * pdf(Binomial(4,b),6-i) for i in 2:3])*a
# batte A all'8vo turno
p72 = sum([pdf(Binomial(4,a),i) * pdf(Binomial(4,b),6-i) for i in 2:4])*a
# batte A al 9no turno
p73 = sum([pdf(Binomial(5,a),i) * pdf(Binomial(4,b),6-i) for i in 2:5])*b
# batte B al 10mo turno
p74 = sum([pdf(Binomial(5,a),i) * pdf(Binomial(5,b),6-i) for i in 1:5])*b
# batte B al 11mo turno
p75 = sum([pdf(Binomial(5,a),i) * pdf(Binomial(6,b),6-i) for i in 0:5])*a
# batte A al 12mo turno
p66 = sum([pdf(Binomial(6,a),i) * pdf(Binomial(6,b),6-i) for i in 0:6])
pda66 = a*b/(1-a-b+2*a*b)
ptot = p70+p71+p72+p73+p74+p75+p66*pda66
@assert 0<=ptot<=1
return ptot
end
aa = range(0.01,0.99, step=0.01)
bb = range(0.01,0.99, step=0.01)
p = surface(aa, bb, ptiebreak,camera=(25,20))
xlabel!("a values")
ylabel!("b values")
contourf(aa, bb, ptiebreak)
xlabel!("a values")
ylabel!("b values")Vediamo ora una delle domande più cruciali nella vita sentimentale: quante persone dovremmo frequentare prima di scegliere quella con cui stabilizzarci per qualcosa di più serio?
È una domanda difficile, ma la matematica può aiutarci trovando una risposta: il 37%. In pratica, di tutte le persone con cui possibilmente potremmo fissare un appuntamento, dovremmo conoscere ma scartare a priori le prime 37%, dalla nostra "scelta definitiva", e poi continuare gli appuntamenti fino ad accasarci con la prima persona che incontriamo e che valutiamo migliore di quelle viste prima; o altrimenti aspettare l'ultima persona, se questa scelta migliore non si presenta.
Perché è una buona strategia?[1] Ovviamente non conviene scegliere la prima persona che incontriamo, per quanto possa sembrare promettente, perché qualcuno di meglio potrebbe apparire più tardi. D'altra parte, non vogliamo nemmeno essere troppo esigenti, perché una volta che rifiutiamo qualcuno molto probabilmente non vorrà tornare da noi qualora ci ripensassimo. Ma quindi perché il 37%? L'idea è solo di massimizzare le probabilità.
Più precisamente: supponiamo di volerci "sistemare" entro un certo periodo. Definiamo \(N\) come il numero di partner che stimiamo di poter incontrare in questo periodo, e che quindi costituiscono i nostri candidati moglie/marito. Supponiamo di poter assegnare a ogni partner un punteggio da 1 a \(N\) per definire il nostro grado di compatibilità[2]. In questo gruppo ci sarà anche la nostra anima gemella, quella giusta, che chiamiamo \(X\) e che avrà quindi il punteggio massimo. Vorremo ora definire una strategia che massimizzi la probabilità di trovarla.
Dato \(N\), quindi, vorremo capire come bilanciare il rischio tra "scegliere un partner iniziale magari non perfetto" e "scartare partner che invece col senno di poi valeva la pena scegliere", e quindi vorremo definire \(M\) come la soglia di persone che delimita i gruppi scarto vs ci penso. La nostra strategia, come sopra anticipato, consiste allora nell'incontrare – ma scartare – le prime \(M\) persone, in modo solo da farci un'idea della qualità delle proposte, e poi scegliere, tra i partner degli appuntamenti successivi, la prima che è migliore della persona migliore incontrata nel gruppo scartato; o altrimenti l'ultima se non si presenta nessuna con questo requisito. Calcolando l'\(M\) ottimale, in funzione di \(N\), risulterà che
\[ M = \frac{N}{e} \approx 0.37 \cdot N \]ovvero \(M\) deve essere circa il 37% di \(N\).
| [1] | "buona strategia" dal punto di vista del risolvere il problema in modo matematico, ma poi naturalmente nelle relazioni umane non ha troppo senso ragionare solo in questi termini razionali. Se per esempio il primo appuntamento ci colpisce non è che dobbiamo rifiutare quel partner solo perché la matematica dice di farlo. |
| [2] | l'assegnazione corretta di questi punteggi esiste ma è a noi ignota, ovviamente, perché il punteggio da assegnare a ogni partner lo sapremo esattamente solo una volta fatti tutti gli appuntamenti. Solo una volta conosciuti tutti potremmo allora dire "lei era la migliore" (punteggio \(N\)), "quest'altra la seconda migliore" (punteggio \(N-1\)), ecc. Durante il percorso invece possiamo solo fare confronti del tipo "lei era peggiore di quella vista prima", "quest'altra era migliore di tutte le altre viste finora", ecc, ovvero confronti relativi. |
using Random
Random.seed!(42)
function M_strategy_score(partners,M,verbose=0)
N = length(partners)
@assert N >= M+1
partners_not_discarded = partners[M+1:end]
best_score_of_discarded = maximum(partners[1:M]) # il punteggio migliore tra le persone del gruppo iniziale
available_scores = max.(partners_not_discarded .- best_score_of_discarded,0)
out = findfirst(u -> u!=0, available_scores)
chosen_index = isnothing(out) ? lastindex(partners_not_discarded) : out
score_out = partners_not_discarded[chosen_index]
found_soul_mate = score_out==N # se score==N abbiamo trovato la nostra anima gemella
if verbose==1
println("We discarded $M partners (from $N partners).")
println("The best score there was $best_score_of_discarded.")
printstyled(partners[1:M]; color=:yellow)
printstyled(partners[M+1:end]; color=:green)
println("\nIn the end we selected partner #$(chosen_index+M) of score $score_out.")
println(found_soul_mate==1 ? "We" : "We didnt have", " found our soulmate.\n")
end
return found_soul_mate # return 0 or 1
end
M = 5
partners = shuffle(1:20)
M_strategy_score(partners,M,1)
# optimal strategy should be discarding M=N/ℯ partners
N = 100
simulations = 30_000
values_tested = ["2","ℯ","3","4","5","6","7"]
M_values = round.(Int64,[N/2 N/ℯ N/3 N/4 N/5 N/6 N/7])
scores = zeros(length(M_values))
for sim in 1:simulations
local partners = shuffle(1:N)
for i in 1:length(M_values)
scores[i] += M_strategy_score(partners,M_values[i])
end
end
scores ./= simulations
println("Testing on N = $N.")
for i in 1:length(M_values)
println("Using M = $(M_values[i]) (that is N/$(values_tested[i])) we got score = $(scores[i])")
end
println("Where score is the probability of having found our soulmate.")
using Plots
plot(bar(scores),xticks=(1:length(values_tested),"N".*"/".*values_tested))We discarded 5 partners (from 20 partners).
The best score there was 19.
[19, 10, 2, 1, 6][12, 9, 8, 11, 15, 7, 4, 17, 20, 16, 5, 3, 14, 18, 13]
In the end we selected partner #14 of score 20.
We found our soulmate.
Testing on N = 100.
Using M = 50 (that is N/2) we got score = 0.3525
Using M = 37 (that is N/ℯ) we got score = 0.37606666666666666
Using M = 33 (that is N/3) we got score = 0.37456666666666666
Using M = 25 (that is N/4) we got score = 0.3532666666666667
Using M = 20 (that is N/5) we got score = 0.3268333333333333
Using M = 17 (that is N/6) we got score = 0.3065333333333333
Using M = 14 (that is N/7) we got score = 0.28196666666666664
Where score is the probability of having found our soulmate.
Problema ispirato da https://plus.maths.org/content/mathematical-dating.
Si consideri una rete di traffico descritta da un grafo \(G = (V, A)\) i cui archi rappresentano le strade di una città e nodi rappresentano gli incroci. I cittadini si spostano ogni giorno da alcune origini verso alcune destinazioni. Sia \(d_k\) la quantità di auto che lasciano punto \(s_k \in V\) e si spostano verso la destinazione \(t_k \in V\). Sia \(\text{FFT}_{ij}\) (free flow time) il tempo di percorrenza della strada (arco) \((i, j)\) se la strada \((i, j)\) è vuota. Il tempo di viaggio su un arco aumenta però naturalmente con il flusso crescente su quell'arco. Una soglia di congestione \(u_{ij}\) viene quindi definita per ogni arco e non deve essere superata dal traffico lungo quella strada.
Modellare il problema di assegnazione del traffico alle strade in modo da ridurre al minimo il tempo totale di viaggio.
using JuMP
using Ipopt
using Plots
using Random
Random.seed!(42)
# Example data for nodes (V), arcs (A), and demands (K)
n_nodes = 5
V = 1:n_nodes # Nodes
A = [(i, j) for i in 1:n_nodes for j in 1:n_nodes if i != j] # Arcs
K = [(1, 4), (2, 3), (2, 4), (2, 5)] # Origin-destination pairs
d = 100*ones(length(K)) # Demand for each origin-destination pair
# Free flow travel time for each arc (FFTi) and congestions u_ij
FFTi = Dict{Tuple{Int, Int}, Float64}()
u = Dict{Tuple{Int, Int}, Float64}()
for arc in A
FFTi[arc] = rand(1:10) # Random travel time between 1 and 10
u[arc] = rand(80:120) # Random congestion threshold between 80 and 120
end
# Create the optimization model
model = Model(Ipopt.Optimizer)
# Define variables: x[k,i,j] for each demand pair (k) and arc (i,j)
@variable(model, x[1:length(K), 1:length(A)] >= 0)
# Objective: Minimize total travel time
@objective(model, Min, sum(FFTi[A[j]] * sum(x[k, j] for k in 1:length(K)) for j in 1:length(A)))
# Flow conservation constraints
for k in 1:length(K)
s_k, t_k = K[k]
for i in V
if i == s_k
@constraint(model, sum(x[k, j] for j in 1:length(A) if A[j][1] == i) - sum(x[k, j] for j in 1:length(A) if A[j][2] == i) == d[k])
elseif i == t_k
@constraint(model, sum(x[k, j] for j in 1:length(A) if A[j][1] == i) - sum(x[k, j] for j in 1:length(A) if A[j][2] == i) == -d[k])
else
@constraint(model, sum(x[k, j] for j in 1:length(A) if A[j][1] == i) - sum(x[k, j] for j in 1:length(A) if A[j][2] == i) == 0)
end
end
end
# Congestion threshold constraints
for j in 1:length(A)
@constraint(model, sum(x[k, j] for k in 1:length(K)) <= u[A[j]])
end
# Solve the model
optimize!(model)
# Extract results
x_opt = value.(x)
# Print results
for k in 1:length(K)
println("Traffic flow for demand pair $(K[k]):")
for j in 1:length(A)
println("Arc $(A[j]): $(x_opt[k, j]) cars")
end
end
# Plot the traffic network with flows
theta = 2 * π * (0:n_nodes-1) / n_nodes
x_plot = 1 * cos.(theta)
y_plot = 1 * sin.(theta)
arc_plot = [(x_plot[i], x_plot[j], y_plot[i], y_plot[j]) for (i, j) in A]
# Generate a list of colors based on the number of demand pairs
colors = [RGB(rand(), rand(), rand()) for _ in 1:length(K)]
# Function to slightly shift points for better visualization
function shift_points(x1, y1, x2, y2, shift_factor)
dx = x2 - x1
dy = y2 - y1
dist = sqrt(dx^2 + dy^2)
shift_x = -dy / dist * shift_factor
shift_y = dx / dist * shift_factor
return (x1 + shift_x, y1 + shift_y, x2 + shift_x, y2 + shift_y)
end
plot()
for j in 1:length(A)
x1, x2, y1, y2 = arc_plot[j]
for k in 1:length(K)
flow = x_opt[k, j]
if flow > 0
sx1, sy1, sx2, sy2 = shift_points(x1, y1, x2, y2, 0.04 * k)
plot!([sx1, sx2], [sy1, sy2], arrow=:arrow, label="",
color=colors[k], linewidth=2,# + flow/50,
alpha=0.6)
mid_x = (sx1 + sx2) / 2
mid_y = (sy1 + sy2) / 2
annotate!(mid_x, mid_y, text("$(round(Int64,flow))", 8, colors[k]))
end
end
end
scatter!(x_plot, y_plot, label="Nodes", legend=:outertopleft)
annotate!(x_plot, y_plot, string.(V))
for k in 1:length(K)
label = "Path $(K[k][1])→$(K[k][2])"
plot!([], [], color=colors[k], label=label)
end
title!("Traffic Network with Flows")
This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.
Number of nonzeros in equality constraint Jacobian...: 160
Number of nonzeros in inequality constraint Jacobian.: 80
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 80
variables with only lower bounds: 80
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 20
Total number of inequality constraints...............: 20
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 20
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.9999960e+00 1.00e+02 7.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 4.1623134e+00 1.00e+02 7.00e+00 -1.0 5.03e+01 - 4.96e-04 4.96e-04h 1
2 2.0988917e+02 9.15e+01 9.43e+00 -1.0 1.00e+02 - 2.23e-04 8.43e-02h 1
3 5.0421487e+02 7.92e+01 6.43e+00 -1.0 9.28e+01 - 3.08e-01 1.34e-01h 1
4 1.9743373e+03 1.75e+01 7.49e+00 -1.0 7.93e+01 - 1.21e-02 7.79e-01h 1
5 1.9924851e+03 1.67e+01 7.04e+00 -1.0 2.70e+01 - 6.14e-02 4.59e-02h 1
6 2.0759046e+03 1.33e+01 7.24e+00 -1.0 2.78e+01 - 2.05e-02 2.03e-01h 1
7 2.1551123e+03 1.03e+01 7.24e+00 -1.0 3.57e+01 - 4.91e-02 2.30e-01h 1
8 2.1629523e+03 9.96e+00 2.80e+00 -1.0 3.03e+01 - 5.85e-01 2.97e-02h 1
9 2.1640031e+03 9.91e+00 3.09e+00 -1.0 9.93e+01 - 2.10e-02 4.79e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.2568808e+03 6.01e+00 4.20e+00 -1.0 1.01e+02 - 1.71e-02 3.93e-01h 1
11 2.3867583e+03 2.02e+00 3.97e+00 -1.0 5.10e+01 - 1.20e-01 6.64e-01h 1
12 2.4369340e+03 7.93e-01 9.18e-01 -1.0 6.10e+00 - 7.83e-01 6.08e-01h 1
13 2.4696395e+03 1.21e-09 3.09e-01 -1.0 4.05e+00 - 5.24e-01 1.00e+00h 1
14 2.4696125e+03 3.22e-10 4.90e-03 -1.0 4.16e+00 - 9.87e-01 1.00e+00f 1
15 2.4643106e+03 1.12e-10 2.00e-07 -1.7 3.21e-01 - 1.00e+00 1.00e+00f 1
16 2.4630097e+03 1.03e-11 4.28e-06 -3.8 6.23e-02 - 1.00e+00 1.00e+00f 1
17 2.4630001e+03 2.84e-14 1.85e-11 -5.7 1.05e-02 - 1.00e+00 1.00e+00f 1
18 2.4630000e+03 1.42e-14 2.66e-14 -8.6 1.51e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 18
(scaled) (unscaled)
Objective...............: 2.4629999790753836e+03 2.4629999790753836e+03
Dual infeasibility......: 2.6645352591003757e-14 2.6645352591003757e-14
Constraint violation....: 1.4210854715202004e-14 1.4210854715202004e-14
Variable bound violation: 9.8526023534162266e-09 9.8526023534162266e-09
Complementarity.........: 3.4143293206596104e-09 3.4143293206596104e-09
Overall NLP error.......: 3.4143293206596104e-09 3.4143293206596104e-09
Number of objective function evaluations = 19
Number of objective gradient evaluations = 19
Number of equality constraint evaluations = 19
Number of inequality constraint evaluations = 19
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 0.021
EXIT: Optimal Solution Found.Problema tratto dagli esercizi del corso di Optimization al polimi
Dato un numero intero positivo \(n\), trovare il poligono con \(n\) lati con diametro (ovvero la distanza massima tra due vertici) minore o uguale a 1, e di area massima.
Un suo esempio di applicazione reale sarebbe massimizzare l'area delimitata da un recinto, imponendo che questo possa essere chiuso in sottorecinti a piacere se necessario, usando staccionate unitarie (ovvero la distanza tra ogni coppia di vertice deve essere ≤1).
using JuMP
using Ipopt
using Plots
# Define the number of vertices
n = 5 # You can change this to the desired number of vertices
# Create the optimization model
model = Model(Ipopt.Optimizer)
# Define variables: r and phi for each vertex
@variable(model, r[1:n] >= 0)
@variable(model, 0 <= phi[1:n] <= 2π)
# Define the objective function (maximize area)
@objective(model, Max,
0.5 * sum(r[i] * r[i+1] * sin(phi[i+1] - phi[i]) for i in 1:n-1) +
0.5 * r[n] * r[1] * sin(2π - phi[n] + phi[1])
)
# Define the constraints for the distance between each pair of vertices
@constraint(model, [i=1:n, j=1:n; j > i], r[i]^2 + r[j]^2 - 2 * r[i] * r[j] * cos(phi[i] - phi[j]) <= 1)
# Define the constraints for ordering angles
@constraint(model, [i=1:n-1], phi[i] <= phi[i+1])
# Solve the model
optimize!(model)
# Get the optimal values for r and phi
r_opt = value.(r)
phi_opt = value.(phi)
# Convert polar coordinates to Cartesian coordinates for plotting
x = r_opt .* cos.(phi_opt)
y = r_opt .* sin.(phi_opt)
# Add the first vertex to close the polygon
x = vcat(x, x[1])
y = vcat(y, y[1])
println("Case n=$n")
for i in 1:length(x)
println("P_$i = (",round(x[i],digits=4), ", ", round(y[i],digits=4),")")
end
# Plot the resulting polygon
plot(x, y, seriestype=:shape, label="Optimal Polygon", fillalpha=0.3, aspect_ratio=:equal, legend=false)
scatter!(x, y, label="Vertices")
title!("Optimal Polygon with n=$n Vertices")
# Function to generate points on a circle
function circle_points(center_x, center_y, radius, num_points=100)
theta = range(0, 2*pi, length=num_points)
cx = center_x .+ radius .* cos.(theta)
cy = center_y .+ radius .* sin.(theta)
return cx, cy
end
# Plot a circle of radius 1 from each vertex
for i in 1:n
circle_x, circle_y = circle_points(x[i], y[i], 1)
plot!(circle_x, circle_y, label=false, linestyle=:solid, color=:red, alpha=0.3)
end
xlims!(minimum(x)-0.2, maximum(x)+0.2)
ylims!(minimum(y)-0.2, maximum(y)+0.2)
Problema tratto dagli esercizi del corso di Optimization al polimi
In quanti modi si può salire una scala lunga \(n\) gradini, potendo salire al massimo 2 scalini per volta? Chiamando \(S(n)\) questo numero di modi in funzione di \(n\), risulta che la soluzione è data da
\[ S(n) = S(n-1) + S(n-2) \]perché moralmente come primo passo dobbiamo sempre fare un passo lungo 1 o lungo 2, indipendentemente da quanto la scala sia alta, e da quei passi poi possiamo osservare di fronte a noi una scala alta \(n-1\) ed \(n-2\) scalini, rispettivamente, e quindi riciclare i calcoli già fatti prima per le altre scale.
Mentre generalizzando il calcolo, supponendo cioè di avere un altro parametro \(k\) che regola quanti scalini risuciamo a salire al massimo con un solo passo, la soluzione diventa
\[ S(n) = S(n-1) + S(n-2) + \ldots + S(n-k) \]dove però ora i vari valori \(S(n-i)\) si devono trovare in modi a volte meno ovvi. Per capire tutto, insieme ad una semplice idea per convertire il calcolo in codice eseguibile da un computer vi rimando naturalmente al video.
using Plots
function S(n,k)
n==0 && return 0
n==1 && return 1
if n<k
return S(n,n)
elseif n==k
return S(n,k-1)+1
else # n>k
return sum(S(n-i,k) for i in 1:k)
end
end
scale = 1:14
modi_1 = S.(scale,1)
modi_2 = S.(scale,2)
modi_3 = S.(scale,3)
modi_4 = S.(scale,4)
modi_5 = S.(scale,5)
plot(scale, modi_1, marker=:circle, line=:solid, label="k=1", legend=:left)
plot!(scale, modi_2, marker=:circle, line=:solid, label="k=2")
plot!(scale, modi_3, marker=:circle, line=:solid, label="k=3")
plot!(scale, modi_4, marker=:circle, line=:solid, label="k=4")
plot!(scale, modi_5, marker=:circle, line=:solid, label="k=5")
plot!(xticks=collect(scale))
xlabel!("Gradini")
ylabel!("Modi")
title!("Comparison of S(n, k) for different k values")Problema ispirato da https://plus.maths.org/content/its-long-way-top.
Due amiche, Astrid e Beatrice, hanno organizzato un picnic. Per pigrizia hanno però deciso di cucinare un'unica pietanza, il riso freddo, e si sono anche accordate che ognuna preparasse quante porzioni riuscisse e che poi le avrebbero condivise ed eventualmente sistemato i conti alla fine.
Astrid ha portato con sé 3 porzioni di riso freddo, mentre Beatrice 5. Sono sul punto di pranzare quando'ecco che appare anche Chiara!, un'altra loro amica, che chiede di unirsi a loro per il picnic, pur non avendo portato niente, ma ovviamente promettendo di saldare il conto alla fine. Astrid e Beatrice accolgono quindi volentieri la nuova arrivata. A questo punto, dividono ogni porzione in tre parti uguali (dove supponiamo che ogni parte valga 1€[3]), generando quindi 9 porzioncine dalle tre porzioni originali di Astrid e 15 dalle cinque di Beatrice, e se le distribuiscono tra loro. Le tre amiche sono molto affamate, e quindi ognuna finisce interamente il riso che le è stato assegnato.
"Gran bella mangiata!" riassume alla fine Chiara, "ma... quanto vi devo? Vediamo: io ho mangiato 8 porzioncine, dalle vostre 3+5 porzioni originali, per quindi un costo di 8€. Di tutto il riso, tu Astrid ne avevi portato i 3/8, mentre Beatrice 5/8, quindi... dovrei dare 3€ ad Astrid e 5€ a Bea, giusto?". "Mi sembra giusto" conferma Astrid. "No invece!" esclama Beatrice, scioccata dall'incompetenza matematica delle sue amiche.
In effetti Beatrice ha ragione. Ma perché? E quale sarebbe la divisione più corretta?
Il problema può essere esteso in forma generale avendo un picnic con \(k\) amiche iniziali, in cui ognuna porta \(n_1, \ldots, n_k\) porzioni di riso, che sono poi raggiunte da altre \(m\) amiche che non hanno portato nulla, ma che alla fine devono pagare una certa somma alle altre amiche.
Problema ispirato da https://plus.maths.org/content/sharing-cakes.
| [3] | questa è solo una semplificazione. Si potrebbe fissare, in modo più verosimile, il costo della porzione inziale (intera), e non della sotto-porzione (divisa), ma questo complicherebbe solo i calcoli. Comunque una volta trovata la soluzione si può convertirla in valori "relativi" e quindi si può ripristinare il valore esatto in base al costo vero della porzione di partenza. |
Tre signori A, B e C hanno fissato un duello a tre armi – un triello quindi – per risolvere una disputa. Erano amici da molto tempo, solo quest'ultima lite ha irrimediabilmente rotto ogni rapporto, quindi tutti conoscono le caratteristiche di tutti: si sa che l'accuratezza, ovvero la probabilità di colpire il bersaglio, di A è 0.3, di C è 0.5, mentre B è infallibile, lo centra sempre.
È stato fissato l'ordine di inizio: A, B, poi C. Ognuno, al proprio turno, deve sparare un colpo decidendo come vuole il suo bersaglio, valutando per esempio quale sia la sua strategia ottimale per sopravvivere. Ovviamente, uno viene ucciso se colpito, e quindi non può eseguire il proprio turno di sparo. Il triello finisce quando un solo uomo rimane vivo.
Se voi foste il signor A, quale strategia utilizzereste?
Immaginate di essere un allevatore che, incrociando varie razze di gallina, è riuscito ad ottenerne una che depone delle uova particolarmente forti. Volete quindi testare di preciso quanto siano resistenti queste uova, e per farlo decidete di lanciarle dai vari piani di un edificio di 100 piani con l'obiettivo di trovare il piano più alto da cui un uovo può cadere senza rompersi.
Se aveste un solo uovo sapreste subito come fare: prima lanciate l'uovo dal primo piano, se non si rompe passate al secondo, poi al terzo, e così via. Ci vorrebbero, nel caso peggiore, 100 lanci per trovare la soluzione.
Ma ora supponiamo che voi abbiate 2 uova. Quale strategia utilizzereste ora per trovare la risposta in modo da minimizzare il caso peggiore del numero di lanci da effettuare?
Supponete infine di avere \(k\) uova. Quale strategia utilizzereste, in questo caso (sempre per minimizzare il caso peggiore del numero di lanci da effettuare)?
Problema ispirato da https://plus.maths.org/content/dropping-eggs.
La funzione che segue (che è stata complicatuccia da scrivere) simula l'esecuzione del lancio di \(k\) uova, usando un jump_size diverso in funzione del numero di uova rimaste. Per esempio, strat_vector = [14,1] corrisponde al caso in cui abbiamo 2 uova, e finché ne abbiamo due saltiamo di 14, poi 13, poi 12 piani, mentre quando ci rimane un solo uovo allora facciamo passi alti sempre 1. Mentre strat_vector = [18,7,1] corrisponde alla simulazione in cui abbiamo 3 uova, e se ne finché ne abbiamo tre effettuiamo salti alti 18, poi 17, 16, ecc; quando ne abbbiamo due salti alti 7, poi 6, 5, ecc; e infine quando ne abbbiamo solo uno salti alti 1.
function the_egg_breaks(testing_floor, critical_floor)
return testing_floor>critical_floor
end
function print_summary(eggs::Integer, floor::Integer, steps::Integer, jump_size::Integer, solution)
println("jump size = $(rpad(jump_size,2)), ",
"🏯=$(rpad(string(floor),3)), ",
"🥚=$(rpad(string(eggs),2)), ",
"🪜=$(rpad(string(steps),2))")
end
function foundable(solution::Vector, verbose=false)
if verbose println(join(map(string, solution))); end
# 2 was the default value
# 1 is for floors where the egg breaks
# 0 is for floors where the egg does not break
return occursin("01",join(map(string, solution))) ||
join(map(string, solution)) == repeat('0',100)
end
function stratk(divs_strat::Vector, critical_floor::Integer, verbose::Bool)
solution = 2 .* ones(Int,100)
eggs = length(divs_strat)
steps = 0
strats = reverse(copy(divs_strat))
previous_floor = 0
current_floor = strats[end]
jump_size = strats[eggs]
while !foundable(solution,verbose) && steps<=100
interval_focus = findlast(x->x=='2',join(map(string, solution))) -
findfirst(x->x=='2',join(map(string, solution))) + 1
eggs<1 && return NaN
jump_size<1 && return NaN
current_floor<1 && return NaN
if isnothing(findlast(x->x=='1',join(map(string,solution))))
max_index = 100
else
max_index = findlast(x->x=='1',join(map(string,solution)))-1
end
current_floor=min(current_floor,max_index)
if the_egg_breaks(current_floor, critical_floor)
eggs -= 1
steps += 1
solution[current_floor:end] .= 1
if verbose print_summary(eggs,current_floor,steps,jump_size,solution); end
current_floor = previous_floor +1
jump_size = eggs>=1 ? strats[eggs] : 1
else
solution[1:current_floor] .= 0
steps += 1
if verbose print_summary(eggs,current_floor,steps,jump_size,solution); end
if eggs>1 jump_size -= 1 end
previous_floor = current_floor
current_floor += jump_size
end
end
return steps
endjulia> stratk([14,1],58,true)
jump size = 14, 🏯=14 , 🥚=2 , 🪜=1
jump size = 13, 🏯=27 , 🥚=2 , 🪜=2
jump size = 12, 🏯=39 , 🥚=2 , 🪜=3
jump size = 11, 🏯=50 , 🥚=2 , 🪜=4
jump size = 10, 🏯=60 , 🥚=1 , 🪜=5
jump size = 1 , 🏯=51 , 🥚=1 , 🪜=6
jump size = 1 , 🏯=52 , 🥚=1 , 🪜=7
jump size = 1 , 🏯=53 , 🥚=1 , 🪜=8
jump size = 1 , 🏯=54 , 🥚=1 , 🪜=9
jump size = 1 , 🏯=55 , 🥚=1 , 🪜=10
jump size = 1 , 🏯=56 , 🥚=1 , 🪜=11
jump size = 1 , 🏯=57 , 🥚=1 , 🪜=12
jump size = 1 , 🏯=58 , 🥚=1 , 🪜=13
jump size = 1 , 🏯=59 , 🥚=0 , 🪜=14
14
julia> stratk([18,7,1],58,true)
jump size = 18, 🏯=18 , 🥚=3 , 🪜=1
jump size = 17, 🏯=35 , 🥚=3 , 🪜=2
jump size = 16, 🏯=51 , 🥚=3 , 🪜=3
jump size = 15, 🏯=66 , 🥚=2 , 🪜=4
jump size = 7 , 🏯=52 , 🥚=2 , 🪜=5
jump size = 6 , 🏯=58 , 🥚=2 , 🪜=6
jump size = 5 , 🏯=63 , 🥚=1 , 🪜=7
jump size = 1 , 🏯=59 , 🥚=0 , 🪜=8
8Babbo Natale ha ricevuto \(n\) letterine dai bambini di tutto il mondo, ciascuna con la lista di regali desiderati. Dopo aver letto attentamente tutte le letterine, e organizzato le spedizioni dei regali, Babbo Natale decide di rispondere personalmente a ciascun bambino, scrivendo un'altra lettera come risposta.
Tuttavia, nella fretta di preparare tutto per le feste, Babbo Natale e i suoi elfi responsabili della logistica si confondono e così, quando sono sul punto di spedire tutto, finiscono per mettere le lettere di risposta nelle varie buste in modo completamente casuale.
Qual è la probabilità che nessun bambino riceva la lettera di risposta che era destinata a lui?
