restart;

with(Student:-Calculus1):

infolevel[Student[Calculus1]] := 1;  # Displays extra details about the calculation

expr := (1 - x)^2 * Sum(k*x^k/(1 - x^k), k = 1 .. infinity);

# Define a rule to use the asymptotic approximation
Rule[rewrite, Sum(k*x^k/(1 - x^k), k = 1 .. infinity) = Pi^2/(6*(1-x)^2)]
    (Limit(expr, x = 1, left));

# Evaluate numerically to check the value
evalf(%);

Creating problem #1
 

restart;
with(plots):

# Definieer de parameters (je kunt deze aanpassen)
a := 1: b := 1: beta := 1: alpha := 1: delta := 1: mu := 1:
t := 0:  # Voor een statische momentopname

# Definieer de functie u in functie-notatie
u := (x, y, z, a, b, beta, alpha, delta, mu) -> (
    (((-1/4 + (((x + y + y + z)*delta)/2 + beta*a + alpha)*mu) *
    sqrt(2*sqrt(1 + 16*((a^2 + b^2)*beta^2 + 2*a*(delta*(y+z) + alpha)*beta + (delta*(y+z) + alpha)^2)*mu^2 +
    8*(-a*beta - delta*(y+z) - alpha)*mu) + 2 + (-8*a*beta - 8*delta*(y+z) - 8*alpha)*mu) *
    sqrt(2*sqrt(1 + 16*((a^2 + b^2)*beta^2 + 2*a*(delta*(x+y) + alpha)*beta + (delta*(x+y) + alpha)^2)*mu^2 +
    8*(-a*beta - delta*(x+y) - alpha)*mu) + 2 + (-8*a*beta - 8*delta*(x+y) - 8*alpha)*mu)) / 8)
);

# Bereken en toon de onafhankelijke variabelen en parameters in u
indets(u);

# Definieer de modulus van u
u_modulus := (x, y, z) -> sqrt(Re(u(x, y, z, a, b, beta, alpha, delta, mu))^2 + Im(u(x, y, z, a, b, beta, alpha, delta, mu))^2);

# 3D-plot van de modulus van u
plot3d(u_modulus(x, y, 0), x=-2..2, y=-2..2, color=blue, axes=boxed, labels=["x", "y", "|u|"], title="Modulus van de M-lump soliton");

restart;

with(plots):

bifdiagram := proc(a1, a2, N::posint, color_choice)
    local bifdiag, a, da, x, pts, cond1, cond2, j;
    
    # Print the meaning of the input parameters
    printf("Bifurcation Diagram Parameters:\n");
    printf("  a1 = %.2f  : Lower bound of parameter 'a'\n", a1);
    printf("  a2 = %.2f  : Upper bound of parameter 'a'\n", a2);
    printf("  N  = %d    : Number of steps between a1 and a2\n", N);
    printf("  Color = %s : Color of the plotted points\n\n", color_choice);
    
    # Define plot conditions with user-defined color
    cond1 := style = point, symbol = point, symbolsize = 8, view = [a1 .. a2, 0 .. 1], color = color_choice;
    cond2 := labels = ["a", "x"], labelfont = [TIMES, ITALIC, 18], axes = boxed;
    
    # Calculate step size
    da := evalf((a2 - a1)/N);
    bifdiag := NULL;
    
    # Loop through parameter values
    for a from a1 by da to a2 do
        x := 0.6;
        
        # Iterate 100 times to reach a stable state
        for j to 100 do
            x := x * a * (1 - x);
        end do;
        
        pts := NULL;
        
        # Collect 100 points for the bifurcation diagram
        for j to 100 do
            x := x * a * (1 - x);
            pts := pts, [a, x];
        end do;
        
        # Append current plot to bifurcation diagram
        bifdiag := bifdiag, plot([pts], cond1, cond2);
    end do;
    
    # Display the final bifurcation diagram
    display([bifdiag]);
end proc:

# Call the function with a custom color (e.g., "red")
bifdiagram(1, 4, 100, red);

Bifurcation Diagram Parameters:
  a1 = 1.00  : Lower bound of parameter 'a'
  a2 = 4.00  : Upper bound of parameter 'a'
  N  = 100    : Number of steps between a1 and a2
  Color = red : Color of the plotted points
 

setup

restart;
with(plots):
with(plottools):

# ✅ Parameters
GammaVal := 1:
omega := 2*Pi:
t0 := 0.7:
tmin := 0:
tmax := 2:
num_frames := 50:

# ✅ Definieer de correcte functies voor f_real en f_imag
f_real := T -> evalf(exp(-GammaVal*(T - t0)^2) * cos(omega*T)):
f_imag := T -> evalf(exp(-GammaVal*(T - t0)^2) * sin(omega*T)):

# ✅ Dynamisch assenbereik bepalen
min_real := evalf(min([seq(f_real(T), T = tmin .. tmax, 0.05)])):
max_real := evalf(max([seq(f_real(T), T = tmin .. tmax, 0.05)])):
min_imag := evalf(min([seq(f_imag(T), T = tmin .. tmax, 0.05)])):
max_imag := evalf(max([seq(f_imag(T), T = tmin .. tmax, 0.05)]));

# ✅ Correcte plaatsing van de projectievlakken
x_offset := evalf(tmax + 1):  
z_offset := evalf(min_imag - 1):  
y_offset := evalf(max_real + 1):  

# ✅ **Toevoegen van de 3D-curve als continue lijn**
space_curve := spacecurve(
    [T, f_real(T), f_imag(T)],
    T = tmin .. tmax,
    color = blue,
    thickness = 3
):

# ✅ **Toevoegen van de projectievlakken**
real_plane := polygon([[tmin, min_real - 1, z_offset], [tmax, min_real - 1, z_offset],
                       [tmax, max_real + 1, z_offset], [tmin, max_real + 1, z_offset]],
                       transparency=0.7, color=lightgray):

imag_plane := polygon([[x_offset, min_real - 1, min_imag], [x_offset, max_real + 1, min_imag],
                       [x_offset, max_real + 1, max_imag], [x_offset, min_real - 1, max_imag]],
                       transparency=0.7, color=lightgray):

time_plane := polygon([[tmin, y_offset, min_imag], [tmax, y_offset, min_imag],
                       [tmax, y_offset, max_imag], [tmin, y_offset, max_imag]],
                       transparency=0.7, color=lightgray):

# ✅ **Animatie van bewegend punt op de 3D-curve**
anim_curve := animate(
    pointplot3d,
    [[[frame, f_real(frame), f_imag(frame)]], color=blue, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het reëel-z vlak**
anim_proj1 := animate(
    pointplot3d,
    [[[frame, f_real(frame), z_offset]], color=red, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het x-imag vlak**
anim_proj2 := animate(
    pointplot3d,
    [[[x_offset, f_real(frame), f_imag(frame)]], color=green, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het y-imag vlak**
anim_proj3 := animate(
    pointplot3d,
    [[[frame, y_offset, f_imag(frame)]], color=purple, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Combineer alle elementen (curve, projectievlakken en animaties)**
animation := display(
    [space_curve, real_plane, imag_plane, time_plane, anim_curve, anim_proj1, anim_proj2, anim_proj3],
    view = [tmin .. x_offset, min_real .. y_offset, z_offset .. max_imag],
    axes = normal,
    orientation = [45, 45]
):

# ✅ Toon de animatie
animation

restart:
# Definieer symbolische variabelen
assume(a, real); assume(b, real); assume(alpha, real); assume(beta, real);
assume(x, real); assume(y, real); assume(t, real);

# Stap 1: Schrijf theta expliciet uit
theta := x + a*y - (alpha + beta*a/(a^2+b^2))*t + I*(b*y - beta*b*t/(a^2+b^2));

# Stap 2: Scheid reële en imaginaire delen
Re_theta := simplify(Re(theta));
Im_theta := simplify(Im(theta));

# Stap 3: Stel Re(theta)=0 en Im(theta)=0 op piekpositie
eq1 := Re_theta = 0;
eq2 := Im_theta = 0;

# Stap 4: Los vergelijking eq2 op naar t
t_sol := solve(eq2, t);

# Stap 5: Vervang t in vergelijking eq1 om expliciete baanvergelijking te vinden
baan_eq := simplify(subs(t = t_sol, eq1));

# Oplossen naar y geeft de expliciete baanvergelijking
baan := solve(baan_eq, y);

Warning, solve may be ignoring assumptions on the input variables.

Warning, solve may be ignoring assumptions on the input variables.

restart;

f := proc(N)
    uses combinat;
    local M, coeff, all_terms, chosen_subsets, s, matchings, matching, pair, sorted_pair, B_product, theta_indices, theta_product, term, i, k;
    all_terms := 0;
    for M from 0 to floor(N/2) do
        ## coeff := 1/(M! * 2^M);  
        coeff := 1;
        if M = 0 then
            all_terms := all_terms + coeff * mul(theta[i], i=1..N);
        else
            chosen_subsets := choose([$1..N], 2*M);
            for s in chosen_subsets do
                matchings := setpartition(s, 2);
                for matching in matchings do
                    B_product := 1;
                    for pair in matching do
                        sorted_pair := sort(pair);
                        B_product := B_product * B[parse(cat(op(sorted_pair)))];
                    end do;
                    theta_indices := {$1..N} minus {op(s)};
                    theta_product := mul(theta[k], k = theta_indices);
                    term := coeff * B_product * theta_product;
                    all_terms := all_terms + term;
                end do;
            end do;
        end if;
    end do;
    return expand(all_terms);
end proc:

F__2 = f(2);

F__4 = f(4)

# Voorbeeld voor N=6
F__6 = f(6);