| > |
# Translated and corrected version: one plot that was not visible is now fixed
restart;
with(PDEtools):
with(Physics):
with(DEtools):
with(plots):
Setup(mathematicalnotation = true):
declare(y(x), prime=x);
# ------------------------------------------------------------
# 1. Helper functions (previously used)
# ------------------------------------------------------------
detect_order := proc(eq, y, x)
local derivs;
derivs := indets(eq, specfunc(diff));
if derivs = {} then return 0; end if;
return max(map(d -> nops([op(d)])-1, derivs));
end proc:
EL_equation := proc(L, y, x)
local dLdy, dLdyprime;
dLdy := diff(L, y(x));
dLdyprime := diff(L, diff(y(x), x));
return simplify(dLdy - diff(dLdyprime, x));
end proc:
detect_type := proc(L)
local partials;
partials := indets(L, 'diff'(anything, anything));
if partials <> {} and ormap(d -> nops([op(d)]) > 2, partials) then
return "field theory (PDE)";
elif nops(indets(L, specfunc(diff))) > 1 then
return "multiple functions (system)";
else
return "classical (1 function, ODE)";
end if;
end proc:
is_linear := proc(eq, y, x)
local vars, max_order, k;
max_order := detect_order(eq, y, x);
vars := {y(x)} union {seq(diff(y(x), x$k), k=1..max_order)};
try
if not type(eq, `=`) then
eq := eq;
else
eq := lhs(eq) - rhs(eq);
end if;
return andmap(term -> degree(term, vars) <= 1, [op(eq)]);
catch:
return false;
end try;
end proc:
historical_stage := proc(expr, y, x, source)
local has_first, has_second, has_x, L_simple, L_brach, L_catenoid, L_beam;
L_simple := sqrt(1 + diff(y(x), x)^2);
L_brach := sqrt((1 + diff(y(x), x)^2)/y(x));
L_catenoid := y(x) * sqrt(1 + diff(y(x), x)^2);
L_beam := diff(y(x), x$2)^2;
if source = "L" then
if expr = L_simple then
return "Pre-Euler (Fermat, 1662) – shortest time / shortest path";
elif expr = L_brach then
return "Pre-Euler (Bernoulli, 1696) – brachistochrone / tautochrone";
elif expr = L_catenoid then
return "Euler (1744) – minimal surface of revolution (catenoid)";
elif expr = L_beam then
return "Post-Euler (19th century) – Euler-Bernoulli beam";
end if;
has_first := has(expr, diff(y(x), x));
has_second := has(expr, diff(y(x), x$2));
has_x := has(expr, x) and not has(expr, diff(y(x), x));
if has_second then
return "Post-Euler (19th century) – higher-order Lagrangian";
elif has_first and has_x then
return "Lagrange (1755) – L = L(x, y, y')";
elif has_first then
return "Euler (1744) – L = L(y, y') autonomous";
else
return "Pre-Euler (17th century) – L without y' (static)";
end if;
else # source = "EL"
has_first := has(expr, diff(y(x), x));
has_second := has(expr, diff(y(x), x$2));
has_x := has(expr, x);
if has_second then
return "Post-Euler – second-order (or higher) ODE";
elif has_first and has_x then
return "Lagrange – first-order ODE with explicit x";
elif has_first then
return "Euler – first-order autonomous ODE";
else
return "Pre-Euler – algebraic equation";
end if;
end if;
end proc:
# ------------------------------------------------------------
# 2. Overview table of famous historical examples
# ------------------------------------------------------------
print_famous_examples_overview := proc()
printf("\n");
printf("====================================================================================================================================\n");
printf("| FAMOUS HISTORICAL EXAMPLES FROM THE CALCULUS OF VARIATIONS |\n");
printf("====================================================================================================================================\n");
printf("| Historical example | Lagrangian L | Euler-Lagrange equation | ODE type (odeadvisor) |\n");
printf("|---------------------------------------|---------------------------------------------|---------------------------------------------|--------------------------------|\n");
printf("| Fermat (1662) – shortest time/path | L = sqrt(1+y'^2) | y'' = 0 | _linear, _second_order, _quadrature |\n");
printf("| Bernoulli (1696) – brachistochrone | L = sqrt((1+y'^2)/y) | 2y y'' + (y')^2 + 1 = 0 | _second_order, _missing_x |\n");
printf("| Huygens (1673) – tautochrone | L = sqrt((1+y'^2)/y) (same) | 2y y'' + (y')^2 + 1 = 0 | _second_order, _missing_x |\n");
printf("| Euler (1744) – catenoid (min. area) | L = y * sqrt(1+y'^2) | y y'' - (y')^2 - 1 = 0 | _second_order, _reducible, _missing_x |\n");
printf("| Lagrange (1755) – explicit x | L = e^x * (y')^2 | y'' + y' = 0 | _linear, _second_order, _exact |\n");
printf("| Geodesic (Lagrange, 1760) – shortest path| L = sqrt(1+y'^2+y^2) | (1+y'^2) y'' - y (1+y'^2) - y (y')^2 = 0? | _second_order, _nonlinear |\n");
printf("| Post-Euler – Euler-Bernoulli beam | L = (y'')^2 | y'''' = 0 | _linear, _fourth_order, _quadrature |\n");
printf("====================================================================================================================================\n");
printf("Note: The tautochrone has the same Lagrangian as the brachistochrone; the difference lies in the boundary conditions.\n");
printf("For the geodesic a simple non-flat metric is chosen: L = sqrt(1 + y'^2 + y^2).\n\n");
end proc:
# ------------------------------------------------------------
# 3. Main analysis procedure (with odeadvisor)
# ------------------------------------------------------------
analyze_Lagrangian := proc(L, y, x,
{ classify_from::string := "L",
show_ode::boolean := false,
show_timeline::boolean := false })
local eq, type_prob, order_prob, lin, history, source_expr, ode_type;
if show_timeline then
print_Lagrangian_timeline();
end if;
type_prob := detect_type(L);
eq := EL_equation(L, y, x);
if show_ode then
printf("Euler-Lagrange equation: %a\n", eq);
end if;
try
ode_type := DEtools:-odeadvisor(eq, y(x));
catch:
ode_type := "cannot determine (not an ODE or too complex)";
end try;
if classify_from = "L" then
source_expr := L;
else
source_expr := eq;
end if;
order_prob := detect_order(eq, y, x);
lin := is_linear(eq, y, x);
history := historical_stage(source_expr, y, x, classify_from);
return table([
"Lagrangian" = L,
"Problem type" = type_prob,
"Euler-Lagrange equation" = eq,
"ODE type (odeadvisor)" = ode_type,
"Order of ODE" = order_prob,
"Linear?" = lin,
"Classification based on" = classify_from,
"Historical phase" = history
]);
end proc:
print_Lagrangian_timeline := proc()
printf("\n");
printf("================================================================================\n");
printf("| HISTORICAL TIMELINE OF THE LAGRANGIAN |\n");
printf("================================================================================\n");
printf("| Period | Figures | Characteristic Lagrangian |\n");
printf("|------------------|--------------------------|-------------------------------|\n");
printf("| Pre-Euler (17th) | Fermat, Bernoulli | L = sqrt(1+y'^2) (Fermat) |\n");
printf("| | | L = sqrt((1+y'^2)/y) (Bernoulli) |\n");
printf("| Euler (1744-60) | Leonhard Euler | L = L(y, y') autonomous |\n");
printf("| Lagrange (1755) | J.-L. Lagrange | L = L(x, y, y') |\n");
printf("| Post-Euler (19th)| Hamilton, Jacobi, Noether| L with higher derivatives |\n");
printf("| Modern | Einstein, Feynman | Lagrangian density for fields |\n");
printf("================================================================================\n");
printf("\n Note: The pre-Euler examples of Fermat and Bernoulli are shown here.\n\n");
end proc:
# ------------------------------------------------------------
# 4. START: show overview of famous examples
# ------------------------------------------------------------
print_famous_examples_overview();
# ------------------------------------------------------------
# 5. EXTRA: Isoperimetric Dido problem (max area, fixed arc length)
# ------------------------------------------------------------
printf("\n");
printf("********************************************************************************\n");
printf("* THE ISOPERIMETRIC PROBLEM OF DIDO (9th century BC) *\n");
printf("* Maximize the area under a curve with given arc length *\n");
printf("* and a straight coastline (the x-axis). *\n");
printf("********************************************************************************\n\n");
# Define the functional with Lagrange multiplier lambda
# We want to maximize: A = ∫ y dx subject to ∫ sqrt(1+y'^2) dx = L (fixed)
# The Lagrangian becomes: L = y + λ * sqrt(1 + y'^2)
lambda := 'lambda': # symbol for the multiplier
L_Dido := y(x) + lambda * sqrt(1 + diff(y(x), x)^2);
printf("Lagrangian with Lagrange multiplier λ: L = y + λ * sqrt(1 + y'^2)\n");
# Euler-Lagrange equation
EL_Dido := EL_equation(L_Dido, y, x);
printf("Euler-Lagrange equation: %a\n", EL_Dido);
# Simplify the equation (it can be reduced to a circle equation)
# Maple can classify the ODE:
try
ode_type_Dido := DEtools:-odeadvisor(EL_Dido, y(x));
printf("ODE type (odeadvisor): %a\n", ode_type_Dido);
catch:
printf("ODE type could not be determined.\n");
end try;
# Solve the ODE (general solution)
sol_general := dsolve(EL_Dido, y(x));
printf("General solution: %a\n", sol_general);
# The solution is a circle: (x - C1)^2 + (y - C2)^2 = λ^2
# We now choose a concrete example: arc length L = 2, coastline from x=-1 to x=1,
# and the curve touches the coast at points (-1,0) and (1,0).
# Then the radius is R = λ, and the centre is (0, R) (circle through (-1,0) and (1,0)).
# Equation: x^2 + (y - R)^2 = R^2 => y = R - sqrt(R^2 - x^2)
# Arc length (half circle) = πR = L = 2 => R = 2/π.
L_length := 2; # given arc length
R := L_length / Pi;
# The optimal curve:
# Correcte definitie van de optimale kromme (halve cirkel op de x-as)
# Correct definition of the optimal curve (semicircle on the x-axis)
y_opt := x -> sqrt(R^2 - x^2);
# Plot of the semicircle
plot_opt := plot(y_opt(x), x=-R..R, color=red, thickness=2, scaling=constrained,
title=sprintf("Optimal circular arc for Dido problem (L = %d, R = %.3f)", L_length, evalf(R)),
labels=["x", "y"], labeldirections=[horizontal, vertical]):
# Coastline (blue dotted line) from -R to R (endpoints exactly touch the arc)
plot_coast := plot(0, x=-R..R, color=blue, thickness=2, linestyle=dot):
# Display
display([plot_opt, plot_coast], axes=boxed);
# Show the computed area
area := int(y_opt(x), x=-R..R);
printf("Enclosed area = %a (exact) ≈ %a\n", area, evalf(area));
# ------------------------------------------------------------
# 6. Analysis of the Dido problem with the existing procedure
# ------------------------------------------------------------
printf("\n--- Analysis via the general procedure (classification) ---\n");
res_Dido := analyze_Lagrangian(L_Dido, y, x, classify_from="L", show_ode=true, show_timeline=false):
print(res_Dido);
# ------------------------------------------------------------
# 7. Extra example: free closed curve (no coastline)
# ------------------------------------------------------------
printf("\n********************************************************************************\n");
printf("* Variant: closed curve (no coastline) – the circle as optimal shape *\n");
printf("********************************************************************************\n");
printf("For a completely closed curve of length L the optimal shape is a circle.\n");
printf("Area = L^2/(4π). For L=2 this gives area = %a.\n", evalf(2^2/(4*Pi)));
printf("This is larger than the area with a coastline (%a), because the coastline restricts freedom.\n", evalf(area));
# ------------------------------------------------------------
# End of code
# ------------------------------------------------------------
printf("\nEnd of demonstration.\n");
|
