P190
projecteuler.net

Maximising a Weighted Product

ℹ️Published on Friday, 18th April 2008, 10:00 pm; Solved by 4586;
Difficulty rating: 50%

Let $S_m = (x_1, x_2, \dots , x_m)$ be the $m$-tuple of positive real numbers with $x_1 + x_2 + \cdots + x_m = m$ for which $P_m = x_1 \cdot x_2^2 \cdot \cdots \cdot x_m^m$ is maximised.

For example, it can be verified that $\lfloor P_{10}\rfloor = 4112$ ($\lfloor \, \rfloor$ is the integer part function).

Find $\sum\limits_{m = 2}^{15} \lfloor P_m \rfloor$.



Soluzione

Esplorazione del problema

Come sempre, per approcciare un problema complicato, possiamo partire in piccolo e osservare cosa accade in casi più semplici; quindi qui per esempio analizzando i casi m=2m=2 e m=3m=3. Giusto perché è carino, e perché Julia permette di farlo molto facilmente, possiamo anche raccontare graficamente questa prima esplorazione.

using Plots
using JSON

Questo è quello che si osserva nel caso di m=2m=2:

nout = 100
x1s = range(0,2,length=nout)
vals = zeros(nout)
for i in 1:nout
	x1 = x1s[i]
	x2 = 2-x1
	vals[i] = x1 * x2^2
end
plot(x1s,vals,label="")
max_idx = argmax(vals)
x1_max = x1s[max_idx]
vline!([x1_max],label="")
println("Max value ($(maximum(vals))) reached at index $max_idx, ie for")
println("x1 = ", x1_max)
println("x2 = ", 2-x1_max)
Max value (1.1851851851851853) reached at index 34, ie for
x1 = 0.6666666666666666
x2 = 1.3333333333333335

vera soluzione:

2-element Vector{Float64}:
 0.6666666666666666
 1.3333333333333333

mentre questo per m=3m=3:

nout = 200
x1s = range(0,3,length=nout)
x2s = range(0,3,length=nout)
vals = zeros(nout,nout)
for i in 1:nout
	x1 = x1s[i]
	for j in 1:nout
		x2 = x2s[j]
		if x1+x2 <= 3
			x3 = 3-x1-x2
			vals[i,j] = x1 * x2^2 * x3^3
		end
	end
end
max_idx = argmax(vals)
x1_max = x1s[max_idx[1]]
x2_max = x1s[max_idx[2]]
plot3d(x1s,x2s,vals',st=:surface, label="")
plot3d!(repeat([x1_max],2),repeat([x2_max],2),[0,maximum(vals)+0.5],label="")
println("Max value ($(maximum(vals))) reached at $max_idx, ie for")
println("x1 = ", x1_max)
println("x2 = ", x2_max)
println("x3 = ", 3-x1_max-x2_max)
Max value (1.6873721655472487) reached at CartesianIndex(34, 67), ie for
x1 = 0.49748743718592964
x2 = 0.9949748743718593
x3 = 1.507537688442211

vera soluzione:

3-element Vector{Float64}:
 0.5
 1.0
 1.5

Questo approccio, come detto, può servire per ambientarci nel problema, e sarebbe ad esempio utile se potessimo già da questi esempi derivare o intuire un qualche pattern rilevante. Tuttavia, non sembra molto il caso dai numeri derivati qui sopra, quindi possiamo passare direttamente alla trattazione matematica e rigorosa del problema.

Soluzione effettiva

L'idea è che stiamo affrontando un problema di ottimizzazione per una funzione multivariata f(x):RnRf({\bf x}) : \mathbb R^n\to \mathbb R, soggetta al vincolo g(x)=0g({\bf x})=0, dove anche g(x):RnRg({\bf x}): \mathbb R^n\to \mathbb R: un contesto che sembra proprio chiamare in causa i moltiplicatori di Lagrange! Quando ci ho pensato non volevo crederci, e soprattutto non volevo rimestare nel torbido per recuperare i miei appunti di Analisi 2. A malincuore ho dovuto però farlo, e in effetti con quelli si arriva tranquillamente (e anche molto velocemente) alla soluzione.

I moltiplicatori di Lagrange permettono infatti di convertire un problema di ottimizzazione (inteso come ricerca di punti di massimo/minimo) di una funzione ff soggetta a un vincolo[1] gg in un altro problema di ottimizzazione, stavolta libero, quindi non soggetto a vincoli, di un'altra funzione L\mathcal{L}, al costo di aggiungere una variabile λ\lambda. Questa nuova funzione sarebbe

L(x,λ)=f(x)+λg(x)\mathcal{L}({\bf x},\lambda) = f({\bf x}) + \lambda g({\bf x})

i cui punti di massimo/minimo sono ora molto facili da individuare: basta infatti usare il classico metodo dell'imporre la "derivata uguale a zero", che in questo contesto multivariato si traduce nel chiedere che il gradiente della funzione sia uguale a zero, ovvero L(x,λ)=0\nabla \mathcal{L}({\bf x},\lambda) = {\bf 0}, per trovare i candidati punti di massimo/minimo. Questa richiesta, vista la definizione di L\mathcal{L}, equivale in realtà a

{f(x)=λg(x)g(x)=0\begin{cases}\nabla f({\bf x}) = -\lambda g({\bf x}) \\ g({\bf x}) = 0\end{cases}

[1] Questo metodo permette anche di gestire più vincoli ovviamente, e in tal caso si introdurrebbe una nuova variabile per vincolo, ma per non appesantire la spiegazione mi sono concentrato sul caso del problema dove c'è un vincolo solo.
Nel nostro caso, il problema è

maxf(x)=i=1mxiis.t.g(x)=i=1mxim=0 \max \, f({\bf x}) = \prod_{i=1}^m {x_i}^i \quad \text{s.t.}\quad g({\bf x}) = \sum_{i=1}^m x_i -m = 0

(s.t.\text{s.t.} sta per "subject to", cioè "soggetto al vincolo") Però, pensando di dover derivare ff per calcolare f(x)\nabla f({\bf x}), questa formulazione non sembra molto comoda, dato che derivando rispetto ad xix_i rimarrebbe comunque un termine che comprende tutte le altre variabili xjx_j. Un'idea è quindi di riformulare il problema applicando a ff un logaritmo: dopotutto, massimizzare f(x)f({\bf x}) è equivalente a massimizzare lnf(x)\ln f({\bf x}), dato che il logaritmo è una funzione crescente. Passando quindi a trattare il problema

maxh(x):=lnf(x)=i=1mixis.t.g(x)=i=1mxim=0 \max \, h({\bf x}):= \ln f({\bf x}) = \sum_{i=1}^m i\cdot{x_i} \quad \text{s.t.}\quad g({\bf x}) = \sum_{i=1}^m x_i -m = 0

e impostando il sistema (2) si ricava

{ixi=λxi=m{xi=iλxi=mλ=m+12xi=2im+1\begin{cases}\frac{i}{x_i} = -\lambda \\ \sum x_i = m \end{cases} \implies \begin{cases}x_i = -\frac{i}{\lambda} \\ \sum x_i = m \end{cases} \implies \lambda = -\frac{m+1}{2} \implies \boxed{x_i = \frac{2i}{m+1}}

da cui segue in poche righe, o con un minimo di impegno anche una sola, la soluzione del problema:

Int.(sum([floor(prod([(2i/(m+1))^i for i in 1:m])) for m in 2:15]))
371048281

Questi sarebbero i valori di xix_i generati per ogni caso:

xs(m) = [(2i/(m+1)) for i in 1:m]
for i in 2:15
	println("$i => ", round.(xs(i),digits=2))
end
2 => [0.67, 1.33]
3 => [0.5, 1.0, 1.5]
4 => [0.4, 0.8, 1.2, 1.6]
5 => [0.33, 0.67, 1.0, 1.33, 1.67]
6 => [0.29, 0.57, 0.86, 1.14, 1.43, 1.71]
7 => [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75]
8 => [0.22, 0.44, 0.67, 0.89, 1.11, 1.33, 1.56, 1.78]
9 => [0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8]
10 => [0.18, 0.36, 0.55, 0.73, 0.91, 1.09, 1.27, 1.45, 1.64, 1.82]
11 => [0.17, 0.33, 0.5, 0.67, 0.83, 1.0, 1.17, 1.33, 1.5, 1.67, 1.83]
12 => [0.15, 0.31, 0.46, 0.62, 0.77, 0.92, 1.08, 1.23, 1.38, 1.54, 1.69, 1.85]
13 => [0.14, 0.29, 0.43, 0.57, 0.71, 0.86, 1.0, 1.14, 1.29, 1.43, 1.57, 1.71, 1.86]
14 => [0.13, 0.27, 0.4, 0.53, 0.67, 0.8, 0.93, 1.07, 1.2, 1.33, 1.47, 1.6, 1.73, 1.87]
15 => [0.12, 0.25, 0.38, 0.5, 0.62, 0.75, 0.88, 1.0, 1.12, 1.25, 1.38, 1.5, 1.62, 1.75, 1.88]

Last modified: May 16, 2025. Website built with Franklin.jl and the lovely Julia programming language.