restart;
with(plots):

# ----------------------------------------------------------------------
# PATHOLOGICAL EXAMPLE: STRICTLY CONVEX QUADRATIC FUNCTIONAL
# ----------------------------------------------------------------------
# Functional:  J[x] = ∫_0^{2π} [ 1/2 (dx/dt)^2 + 1/2 x^2 - x·f(t) ] dt
# Euler-Lagrange:  x'' - x = -f(t)
# Boundary conditions: x(0)=0, x(2π)=0
# Forcing (non‑trivial, with three frequencies):
f := t -> sin(t) + 0.3*sin(2*t) + 0.1*sin(3*t);

# ----------------------------------------------------------------------
# 1. EXACT SOLUTION OF THE EULER-LAGRANGE ODE
# ----------------------------------------------------------------------
ode := diff(x(t), t, t) - x(t) = -f(t);
bc  := x(0)=0, x(2*Pi)=0;
sol_exact := dsolve({ode, bc}, x(t)):
x_exact := unapply(rhs(sol_exact), t):

# ----------------------------------------------------------------------
# 2. RITZ METHOD (DIRECT MINIMIZATION) USING A SINE BASIS
#    Because the functional is quadratic, the minimum satisfies a linear
#    system:  A c = b, where
#      A_{nm} = ∫ (φ_n' φ_m' + φ_n φ_m) dt,   φ_n = sin(n t)
#      b_n   = ∫ φ_n f dt
# ----------------------------------------------------------------------
N := 5;   # number of basis functions
# Define basis functions
phi := n -> t -> sin(n*t);

# Compute matrix A and vector b via numerical integration
A := Matrix(N, N, datatype=float[8]):
b := Vector(N, datatype=float[8]):
for n from 1 to N do
    for m from 1 to N do
        integrand_A := diff(phi(n)(t), t)*diff(phi(m)(t), t) + phi(n)(t)*phi(m)(t);
        A[n,m] := evalf(Int(integrand_A, t=0..2*Pi, method=_d01ajc));
    end do;
    integrand_b := phi(n)(t)*f(t);
    b[n] := evalf(Int(integrand_b, t=0..2*Pi, method=_d01ajc));
end do:

# Solve the linear system
c_opt := LinearAlgebra:-LinearSolve(A, b):
coeff_opt := [seq(c_opt[n], n=1..N)];

# Minimal value of J (optional)
J_min := 1/2 * (c_opt^%T . A . c_opt) - (b^%T . c_opt);

# Construct the Ritz approximation
x_ritz := unapply( add(coeff_opt[n]*sin(n*t), n=1..N), t );

# ----------------------------------------------------------------------
# 3. PLOT: COMPARE EXACT SOLUTION AND RITZ APPROXIMATION
# ----------------------------------------------------------------------
p_exact := plot(x_exact(t), t=0..2*Pi, color=blue, thickness=2, legend="Exact (Euler-Lagrange)");
p_ritz  := plot(x_ritz(t), t=0..2*Pi, color=red, linestyle=dash, thickness=2, legend="Ritz (direct min, N=5)");
display(p_exact, p_ritz,
        title="Equivalence: Euler-Lagrange ODE ↔ Direct minimization",
        labels=["t", "x(t)"],
        size=[800,400]);

# ----------------------------------------------------------------------
# 4. OUTPUT WITH EDUCATIONAL EXPLANATION
# ----------------------------------------------------------------------
printf("============================================================\n");
printf("PATHOLOGICAL EXAMPLE: CONVEX QUADRATIC FUNCTIONAL\n");
printf("============================================================\n");
printf("Functional: J[x] = ∫ (½ x'² + ½ x² - x·f(t)) dt\n");
printf("f(t) = sin(t) + 0.3 sin(2t) + 0.1 sin(3t)\n");
printf("Boundary conditions: x(0)=0, x(2π)=0\n\n");

printf("RESULTS:\n");
printf("Exact solution (Euler-Lagrange) computed.\n");
printf("Ritz coefficients (N=%d):\n", N);
for n from 1 to N do
    printf("  c[%d] = %.8f\n", n, coeff_opt[n]);
end do;
printf("Minimal J = %.8f\n", J_min);
printf("\nThe plot shows that both approaches coincide.\n\n");

printf("============================================================\n");
printf("WHAT IS ‘PATHOLOGICAL’ ABOUT THIS EXAMPLE?\n");
printf("============================================================\n");
printf("1. The functional is STRICTLY CONVEX (quadratic, positive definite).\n");
printf("   According to convex analysis (Rockafellar, Hiriart-Urruty/Lemarechal),\n");
printf("   a convex functional has at most one minimum, and the necessary\n");
printf("   condition (Euler-Lagrange) is also sufficient. This example confirms\n");
printf("   that the ODE solution and the direct minimization coincide exactly.\n\n");
printf("2. The Euler-Lagrange equation is LINEAR (x'' - x = -f(t)).\n");
printf("   In most physical problems, the calculus of variations is nonlinear.\n");
printf("   The linearity makes this example ‘artificial’ and therefore ideal\n");
printf("   for teaching: the exact solution is available in closed form.\n\n");
printf("3. Even though the linear ODE with constant coefficients has exponentially\n");
printf("   growing/decaying solutions, the boundary conditions x(0)=x(2π)=0 and\n");
printf("   the choice of the period 2π force the solution to be PERIODIC.\n");
printf("   This is a surprising effect that becomes clearly visible only in a\n");
printf("   ‘pathological’ example.\n\n");
printf("4. The Ritz method uses only sines (which satisfy the BCs). The exact\n");
printf("   solution also contains cosines and exponential terms, but those are\n");
printf("   suppressed by the boundary conditions. Nevertheless, the Ritz\n");
printf("   approximation with N=5 converges almost exactly to the analytical\n");
printf("   solution. This illustrates the power of direct minimization, even\n");
printf("   with a ‘pathological’ choice of basis.\n\n");
printf("In short: the example is ‘pathological’ because it is UNNATURALLY SIMPLE\n");
printf("(linear, convex, exactly solvable), but precisely that simplicity allows\n");
printf("us to demonstrate clearly the theoretical equivalence between the\n");
printf("Euler-Lagrange equation and the minimization problem.\n");
printf("============================================================\n\n");

printf("LITERATURE REFERENCES (as requested):\n");
printf("  Rockafellar, R.T. (1970). Convex Analysis. Princeton Univ. Press.\n");
printf("  Hiriart-Urruty, J.-B. & Lemaréchal, C. (1993). Convex Analysis and\n");
printf("    Minimization Algorithms. Springer.\n");
printf("  Michlin, S.G. (1964). Variation Problems of Mathematical Physics.\n");
printf("    (English translation) Pergamon Press.\n");
printf("============================================================\n");

============================================================
PATHOLOGICAL EXAMPLE: CONVEX QUADRATIC FUNCTIONAL
============================================================
Functional: J[x] = ∫ (½ x'² + ½ x² - x·f(t)) dt
f(t) = sin(t) + 0.3 sin(2t) + 0.1 sin(3t)
Boundary conditions: x(0)=0, x(2π)=0

RESULTS:
Exact solution (Euler-Lagrange) computed.
Ritz coefficients (N=5):
  c[1] = 0.50000000
  c[2] = 0.06000000
  c[3] = 0.01000000
  c[4] = -0.00000000
  c[5] = -0.00000000
Minimal J = -0.81524329

The plot shows that both approaches coincide.

============================================================
WHAT IS ‘PATHOLOGICAL’ ABOUT THIS EXAMPLE?
============================================================
1. The functional is STRICTLY CONVEX (quadratic, positive definite).
   According to convex analysis (Rockafellar, Hiriart-Urruty/Lemarechal),
   a convex functional has at most one minimum, and the necessary
   condition (Euler-Lagrange) is also sufficient. This example confirms
   that the ODE solution and the direct minimization coincide exactly.

2. The Euler-Lagrange equation is LINEAR (x'' - x = -f(t)).
   In most physical problems, the calculus of variations is nonlinear.
   The linearity makes this example ‘artificial’ and therefore ideal
   for teaching: the exact solution is available in closed form.

3. Even though the linear ODE with constant coefficients has exponentially
   growing/decaying solutions, the boundary conditions x(0)=x(2π)=0 and
   the choice of the period 2π force the solution to be PERIODIC.
   This is a surprising effect that becomes clearly visible only in a
   ‘pathological’ example.

4. The Ritz method uses only sines (which satisfy the BCs). The exact
   solution also contains cosines and exponential terms, but those are
   suppressed by the boundary conditions. Nevertheless, the Ritz
   approximation with N=5 converges almost exactly to the analytical
   solution. This illustrates the power of direct minimization, even
   with a ‘pathological’ choice of basis.

In short: the example is ‘pathological’ because it is UNNATURALLY SIMPLE
(linear, convex, exactly solvable), but precisely that simplicity allows
us to demonstrate clearly the theoretical equivalence between the
Euler-Lagrange equation and the minimization problem.
============================================================

LITERATURE REFERENCES (as requested):
  Rockafellar, R.T. (1970). Convex Analysis. Princeton Univ. Press.
  Hiriart-Urruty, J.-B. & Lemaréchal, C. (1993). Convex Analysis and
    Minimization Algorithms. Springer.
  Michlin, S.G. (1964). Variation Problems of Mathematical Physics.
    (English translation) Pergamon Press.
============================================================

# ============================================================================
# DIDACTIC WORKSHEET: DERIVATION OF THE ROTATION MATRIX ABOUT THE X‑AXIS
#
# This worksheet guides you through:
#   1. 2D rotation in the yz‑plane (point on a circle)
#   2. Extension to 3D: rotation of a vector about the x‑axis
#   3. Rotation of a cube about the x‑axis with a dynamic matrix display
#
# Run each section separately to see the theory and the animations.
# ============================================================================

restart;
with(plots):        # for plots and animations
with(plottools):    # for lines, polygons, spheres
with(LinearAlgebra):# for matrix multiplication

# ----------------------------------------------------------------------------
# PART 1: 2D ROTATION IN THE YZ‑PLANE
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 1: 2D rotation in the yz‑plane");
print("We consider a point P = (y, z) in the yz‑plane.");
print("We want to rotate it counter‑clockwise by an angle θ.");
print("");

print("Step 1 – Polar coordinates:");
print("   y = r·cos(φ),   z = r·sin(φ)");
print("where r is the distance from the origin and φ is the angle with the positive y‑axis.");
print("");

print("Step 2 – Rotation by θ adds θ to the angle:");
print("   y' = r·cos(φ+θ),   z' = r·sin(φ+θ)");
print("");

print("Step 3 – Use the addition formulas:");
print("   cos(φ+θ) = cosφ·cosθ - sinφ·sinθ");
print("   sin(φ+θ) = sinφ·cosθ + cosφ·sinθ");
print("");

print("Step 4 – Multiply by r and substitute y = r·cosφ, z = r·sinφ:");
print("   y' = y·cosθ - z·sinθ");
print("   z' = y·sinθ + z·cosθ");
print("");

print("Step 5 – Write as a matrix multiplication:");
print("   [ y' ]   [ cosθ  -sinθ ] [ y ]");
print("   [ z' ] = [ sinθ   cosθ ] [ z ]");
print("The 2×2 matrix is the 2D rotation matrix R_2D(θ).");
print("============================================================================");
print("");

print("Animation: A point starting at (1,0) rotates through a full circle.");
print("Below the plot you see the numerical matrix R_2D(θ).");
print("");

N := 32;  # number of frames
frames_2d := []:
for i from 0 to N do
    theta := i * (2*Pi/N);
    R2 := Matrix([[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]);
    rotated := convert(R2 . Vector([1,0]), list);
    
    # Graphical elements
    p := pointplot([rotated], symbol=solidcircle, symbolsize=20, color="red");
    l := line([0,0], rotated, color="blue", linestyle=dot);
    circ := plot([cos(t), sin(t), t=0..2*Pi], color="gray", thickness=1);
    y_axis := plot([[-1.5,0],[1.5,0]], color="black", thickness=1);
    z_axis := plot([[0,-1.5],[0,1.5]], color="black", thickness=1);
    y_label := textplot([1.6,0.1, "y"], font=["Times",12]);
    z_label := textplot([0.1,1.6, "z"], font=["Times",12]);
    matrix_text := sprintf("R_2D = [%7.3f  %7.3f]\n        [%7.3f  %7.3f]",
                           cos(theta), -sin(theta), sin(theta), cos(theta));
    mat_plot := textplot([-1.2, -1.3, matrix_text], font=["Courier",10]);
    
    frames_2d := [op(frames_2d),
                  display(p, l, circ, y_axis, z_axis, y_label, z_label, mat_plot,
                          scaling=constrained, view=[-1.5..1.5,-1.5..1.5],
                          title=cat("θ = ", sprintf("%.2f", theta), " rad"))];
end do:
anim_2d := display(frames_2d, insequence=true):
print(anim_2d);
print("Click the play button. The point rotates and the matrix updates.");
print("");

# ----------------------------------------------------------------------------
# PART 2: EXTENSION TO 3D – ROTATION ABOUT THE X‑AXIS
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 2: From 2D to 3D – rotation about the x‑axis");
print("In 3D we have coordinates (x, y, z).");
print("When rotating about the x‑axis:");
print("  - The x‑coordinate remains unchanged (the axis is the x‑axis).");
print("  - The y and z coordinates rotate exactly as in the yz‑plane.");
print("Therefore:");
print("   x' = x");
print("   y' = y·cosθ - z·sinθ");
print("   z' = y·sinθ + z·cosθ");
print("");

print("Step 6 – Write as a 3×3 matrix:");
print("   [ x' ]   [ 1    0     0 ] [ x ]");
print("   [ y' ] = [ 0  cosθ -sinθ ] [ y ]");
print("   [ z' ]   [ 0  sinθ  cosθ ] [ z ]");
print("This is the rotation matrix R_x(θ).");
print("============================================================================");
print("");

print("Animation: A red vector initially at (0,1,0) rotates about the x‑axis.");
print("The vector stays in the yz‑plane; its x‑coordinate remains 0.");
print("");

Rx := theta -> Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]]):
frames_vec := []:
N := 32;
for i from 0 to N do
    theta := i * (2*Pi/N);
    v0 := Vector([0,1,0]);
    v_rot := Rx(theta) . v0;
    tip := convert(v_rot, list);
    # Draw vector as a line + a small sphere at the tip
    line_vec := line([0,0,0], tip, color="red", thickness=3);
    tip_sphere := sphere(tip, 0.08, color="red");
    # Coordinate axes as lines
    x_axis := line([-1.5,0,0], [3,0,0], color="black", thickness=1);
    y_axis := line([0,-1.5,0], [0,3,0], color="black", thickness=1);
    z_axis := line([0,0,-1.5], [0,0,3], color="black", thickness=1);
    axis_text := textplot3d([[3.2,0,0,"x"],[0,3.2,0,"y"],[0,0,3.2,"z"]], font=["Times",12]);
    frames_vec := [op(frames_vec),
                   display(line_vec, tip_sphere, x_axis, y_axis, z_axis, axis_text,
                           scaling=constrained, view=[-1.5..3,-1.5..3,-1.5..3],
                           orientation=[70,60], title=cat("θ = ", sprintf("%.2f", theta), " rad"))];
end do:
anim_vec := display(frames_vec, insequence=true):
print(anim_vec);
print("Click the play button. The red vector rotates in the yz‑plane.");
print("");

# ----------------------------------------------------------------------------
# PART 3: ROTATING A CUBE ABOUT THE X‑AXIS
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 3: Rotating a cube about the x‑axis");
print("We apply the matrix R_x(θ) to every vertex of a cube.");
print("The cube has side length 2 and its centre at the origin.");
print("A dynamic caption below the plot shows the current numerical matrix.");
print("============================================================================");
print("");

# Cube definition
cube_vertices := [
    [-1,-1,-1], [ 1,-1,-1], [ 1, 1,-1], [-1, 1,-1],  # bottom face (z = -1)
    [-1,-1, 1], [ 1,-1, 1], [ 1, 1, 1], [-1, 1, 1]   # top face (z =  1)
];
edges := [
    [1,2],[2,3],[3,4],[4,1],  # bottom
    [5,6],[6,7],[7,8],[8,5],  # top
    [1,5],[2,6],[3,7],[4,8]   # vertical
];

draw_cube := proc(vertices)
    local edge, lines, faces;
    lines := seq( line(vertices[edge[1]], vertices[edge[2]], color="Black", thickness=2), edge in edges );
    faces := [
        polygon([vertices[1],vertices[2],vertices[3],vertices[4]], color="Red", transparency=0.5),
        polygon([vertices[5],vertices[6],vertices[7],vertices[8]], color="Green", transparency=0.5),
        polygon([vertices[1],vertices[2],vertices[6],vertices[5]], color="Blue", transparency=0.5),
        polygon([vertices[2],vertices[3],vertices[7],vertices[6]], color="Yellow", transparency=0.5),
        polygon([vertices[3],vertices[4],vertices[8],vertices[7]], color="Cyan", transparency=0.5),
        polygon([vertices[4],vertices[1],vertices[5],vertices[8]], color="Magenta", transparency=0.5)
    ];
    return display(lines, faces);
end:

frames_cube := []:
N := 32;
for i from 0 to N do
    theta := i * (2*Pi/N);
    rotated_vertices := map(v -> convert(Rx(theta) . Vector(v), list), cube_vertices);
    M := Rx(theta);
    matrix_caption := cat(
        "R_x(θ) = \n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[1,1], M[1,2], M[1,3]), "\n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[2,1], M[2,2], M[2,3]), "\n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[3,1], M[3,2], M[3,3])
    );
    frames_cube := [op(frames_cube),
                    display(draw_cube(rotated_vertices),
                            title = cat("θ = ", sprintf("%.2f", theta), " rad"),
                            orientation = [45,45,0],
                            lightmodel = light4,
                            scaling = constrained,
                            axes = normal,
                            labels = ["x","y","z"],
                            view = [-2.5..2.5, -2.5..2.5, -2.5..2.5],
                            caption = matrix_caption
                           )]:
end do:
anim_cube := display(frames_cube, insequence=true):
print(anim_cube);
print("Click the play button. The cube rotates about the x‑axis.");
print("The caption shows the matrix R_x(θ) with current numerical values.");
print("============================================================================");
print("");

print("END OF WORKSHEET.");
print("You have now seen step by step how the rotation matrix about the x‑axis is derived,");
print("from a 2D point rotation, through a 3D vector, to a full cube.");

restart;
with(Fractals[LSystem]):
with(plots):

# Define the L-system generator
LGenerator := proc(n, Init, Angle0, Param, Rules, Color)
    local Iteration;
    
    # Use try-catch to avoid errors
    try
        if n = 0 then
            Iteration := Init;
        else
            Iteration := Iterate(Init, Rules, n);
        end if;
    catch:
        # If there's an error, use the initial string
        Iteration := Init;
    end try;
    
    # Draw the L-system
    LSystemPlot(
        Iteration,
        Param,
        initialangle = Angle0,
        color = Color,
        scaling = constrained,
        axes = none,
        thickness = 2
    );
end proc:

# First L-system (number 17)
Init17 := "X":
Angle17 := 90:

Param17 := [
  "F" = "draw:1",
  "+" = "turn:-25.7",
  "-" = "turn:25.7",
  "[" = "push:position;push:angle",
  "]" = "pop:position;pop:angle"
]:

Rules17 := [
  "F" = "FF",
  "X" = "F[+X][-X]FX"
]:

# Second L-system (number 18)
Init18 := "X":
Angle18 := 90:

Param18 := [
  "F" = "draw:1",
  "+" = "turn:20",
  "-" = "turn:-20",
  "[" = "push:position;push:angle",
  "]" = "pop:position;pop:angle"
]:

Rules18 := [
  "F" = "FF",
  "X" = "F[+X]F[-X]+X"
]:

# Test both L-systems without animation
print("L-system 17, iteration 4:");
display(LGenerator(4, Init17, Angle17, Param17, Rules17, "Red"));

print("L-system 18, iteration 4:");
display(LGenerator(4, Init18, Angle18, Param18, Rules18, "DarkGreen"));

# Create smoother animations with more frames
print("Smooth animation of L-system 17 (iterations 0-7):");
display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],
        insequence = true,
        frames = 64);  # More frames for smoother animation

print("Smooth animation of L-system 18 (iterations 0-7):");
display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],
        insequence = true,
        frames = 64);  # More frames for smoother animation

# Optional: Show all iterations side by side in a grid
print("All iterations of L-system 17 (0-7):");
display(Matrix(2, 4, [seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)]));

print("All iterations of L-system 18 (0-7):");
display(Matrix(2, 4, [seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)]));

# Even smoother version with interpolation (more iterations shown)
print("Extra smooth animation of L-system 17 (more frames):");
display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],
        insequence = true,
        frames = 128);  # Even more frames for extra smoothness

print("Extra smooth animation of L-system 18 (more frames):");
display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],
        insequence = true,
        frames = 128);

# --- Symbols ---
declare(R(xi), phi(xi)):
#assume(alpha::real, beta::real, v::real):  ?????????????????

# --- Imaginary part Eq. (2.4) ---
ImEq :=
alpha*R(xi)*diff(phi(xi),xi$2)
+ 2*alpha*diff(R(xi),xi)*diff(phi(xi),xi)
- v*diff(R(xi),xi)
+ beta*(2*n+1)*R(xi)^(2*n)*diff(R(xi),xi):

# --- Chirp ansatz (Eq. 2.5) ---
phip  := eta_1*R(xi)^(2*n) + eta_2:
phipp := diff(phip, xi):

# --- Substitute ansatz ---
ImEq2 := subs(
  diff(phi(xi),xi)=phip,
  diff(phi(xi),xi$2)=phipp,
  ImEq
):

# --- Divide out R' ---
ImEq3 := factor(ImEq2/diff(R(xi),xi)):

# --- Helper symbol ---
Z := 'Z':
ImEq4 := subs(R(xi)^(2*n)=Z, ImEq3):

# --- Collect ---
C := collect(ImEq4, Z):

# --- Coefficients ---
C_Z := coeff(C, Z, 1):
C_0 := coeff(C, Z, 0):

# --- Solve ---
eta[1] := solve(C_Z = 0, eta_1);
eta[2] := solve(C_0 = 0, eta_2);

#eta1_sol, eta2_sol;