ChaoticCircle := proc(r::positive, factor_x::positive, factor_y::positive, k::positive, delta::numeric, n::posint)
    local i, j, N, t, theta, x_vals, y_vals, points, plot, freqs, phases, wobblefactor;
    uses plots, plottools, RandomTools;

    N := 500;
    theta := [seq(2*Pi*i/N, i = 0..N)];

    # Generate random frequencies and phases
    freqs := [seq(rand(3..7)(), i=1..n)];
    phases := [seq(rand(0.0..evalf(2*Pi))(), i=1..n)];

    # Function to compute wobble factor at t
    wobblefactor := t -> 1 + delta * add(sin(freqs[j]*t + phases[j]), j=1..n);

    # Coordinates with random contour wobble
    x_vals := [seq(k*factor_x*r*wobblefactor(t)*cos(t), t in theta)];
    y_vals := [seq(k*factor_y*r*wobblefactor(t)*sin(t), t in theta)];

    points := [seq([x_vals[i], y_vals[i]], i=1..N+1)];

    # Plot the chaotic shape
    plot := plottools[curve](points, color=orange, thickness=3):
    plots[display](plot,
      axes = none,
        title=sprintf("Chaotic Circle  �� (n=%d waves)", n));
end proc:

ChaoticCircle(1, 5, 2, 4, 0.07, 5);

restart:

# Step 1: Define parameters
assume(a, real): assume(b, real): assume(lambda > 0):
assume(k1, real): assume(k2, real):
assume(r1, real): assume(r2, real):
assume(s1, real): assume(s2, real):

# Step 2: Condition A_12 = 0 → set numerator equal to 0
numerator := -(k1 - k2)^2 + 3*(s1 - s2)^2 + lambda*(r1 - r2)^2 + (a^4 - 6*a^2*b^2 + b^4):
eq := numerator = 0:

# Step 3: Solve the equation for r2
sol := solve(eq, r2):

# Step 4: Look at the difference r2 - r1 (for ± structure)
expr := sol[1] - r1:
expr_simplified := simplify(expr):

# Step 5: Rewrite the result in the same form as equation (14)
A := (a^4 - 6*a^2*b^2 + b^4):
factorized_expr := simplify(sqrt(3)/(4*lambda) * sqrt(-lambda*(k2 - k1)^2*A + 4*lambda*(r1 + s1 - s2))):

# Step 6: Display the final solution
r2_final := r1 + factorized_expr;
r2_alternate := r1 - factorized_expr;

# Display both ± solutions
[r2_final, r2_alternate];

# Replace expanded terms by original factor form
A := a^4 - 6*a^2*b^2 + b^4:
DeltaK := (k1 - k2)^2:
expr_compact := simplify(sqrt(3)/(4*lambda) * sqrt(-lambda*DeltaK*A + 4*lambda*(r1 + s1 - s2))):

# Show the final clean formula
r2_nice := [r1 + expr_compact, r1 - expr_compact];


A := (a^4 - 6*a^2*b^2 + b^4):
DeltaK := (k1 - k2)^2:
B := 4*(r1 + s1 - s2):
r2_compact := [r1 + sqrt(3)/(4*sqrt(lambda))*sqrt(-DeltaK*A + B),
               r1 - sqrt(3)/(4*sqrt(lambda))*sqrt(-DeltaK*A + B)];

 

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

restart:

# Parameters definition
assume(a, real): assume(b, real): assume(lambda > 0):
assume(k1, real): assume(k2, real):
assume(r1, real): assume(s1, real): assume(s2, real):

# Compact terms as in the paper
DeltaK := (k2 - k1)^2:
A := a^4 - 6*a^2*b^2 + b^4:
B := 4*lambda*(r1 + s1 - s2):

# Root term as it appears in equation (14)
Root := -lambda*DeltaK*A + B:

# Final solution in ± form
r2_nice := [ + sqrt(3)/(4*lambda)*sqrt(Root),
             - sqrt(3)/(4*lambda)*sqrt(Root) ];

 

restart;

# Simplification assuming lambda > 0
assume(lambda > 0):

DeltaK := (k2 - k1)^2:
A := a^4 - 6*a^2*b^2 + b^4:
B := 4*(r1 + s1 - s2):  # No lambda in this term anymore

# Extract lambda from the square root
r2_nice := r1 + sqrt(3)/(4*sqrt(lambda)) * sqrt(-DeltaK*A + B),
            r1 - sqrt(3)/(4*sqrt(lambda)) * sqrt(-DeltaK*A + B);

restart;
with(PDEtools): with(plots):

# --------------------------
# 1. Basisparameters
# --------------------------
k[2] := 1: l[2] := 1: r[2] := 0.5:

# Tijdwaarden die je wil onderzoeken
timeList := [-5, -2, 0, 2, 5]:

# G-functie (geëxtraheerd uit lumpoplossing)
t_expr := (tval) -> exp(-(1/k[2])*(tval*k[2]^4 + tval*k[2]^2 + tval*k[2]*l[2] + tval*k[2]*r[2] - tval*l[2]^2 - x*k[2]^2 - y*k[2]*l[2])) +
                    exp(-1) +
                    exp(-(1/k[2])*(tval*k[2]^4 + tval*k[2]^2 + tval*k[2]*l[2] + tval*k[2]*r[2] - tval*l[2]^2 - x*k[2]^2 - y*k[2]*l[2])):

# --------------------------
# 2. Tijdsafhankelijke plot
# --------------------------
plotsList := []:

for T in timeList do
    fxy := t_expr(T):
    fx := diff(fxy, x): fxx := diff(fxy, x$2):
    u := unapply(2*fxx/fxy - 2*(fx/fxy)^2, x, y):
    P := plot3d(u(x, y), x = -20 .. 20, y = -20 .. 20,
                title = cat("Tijdsontwikkeling bij t = ", T),
                labels = ["x", "y", "u(x,y)"],
                axes = boxed, grid = [60,60], color = gold):
    plotsList := [op(plotsList), P]:
end do:

display(plotsList, insequence = true);

# --------------------------
# 3. Contouroppervlak op t = 0
# --------------------------
fxy0 := t_expr(0):
fx0 := diff(fxy0, x): fxx0 := diff(fxy0, x$2):
u0 := unapply(2*fxx0/fxy0 - 2*(fx0/fxy0)^2, x, y):

contourplot3d(u0(x, y), x = -20 .. 20, y = -20 .. 20, contours = 25,
              title = "Contourplot van Lump bij t = 0",
              labels = ["x", "y", "u(x,y)"]);

restart;
with(PDEtools): with(plots):

# ---------------------------------------------------
# Stap 1: Genereer lump-functie als functie van t, x, y
# ---------------------------------------------------
Lump := proc(k2, l2, r2, tval)
    local G, fxy, fx, fxx, u;
    G := exp(-(1/k2)*(tval*k2^4 + tval*k2^2 + tval*k2*l2 + tval*k2*r2 - tval*l2^2 - x*k2^2 - y*k2*l2)) +
         exp(-1) +
         exp(-(1/k2)*(tval*k2^4 + tval*k2^2 + tval*k2*l2 + tval*k2*r2 - tval*l2^2 - x*k2^2 - y*k2*l2)):
    fx := diff(G, x): fxx := diff(G, x$2):
    u := unapply(2*fxx/G - 2*(fx/G)^2, x, y):
    return u:
end proc:

# ---------------------------------------------------
# Stap 2: Parameterwaarden om te variëren
# ---------------------------------------------------
Kvals := [0.5, 1, 2];     # k[2]: invloed op scherpte
Lvals := [-1, 0, 1];      # l[2]: invloed op y-positie
Rvals := [0, 0.5, 1];     # r[2]: invloed op z-as (indirect hier)

# ---------------------------------------------------
# Stap 3: Genereer alle plots automatisch
# ---------------------------------------------------
plotsList := []:
for k2 in Kvals do
  for l2 in Lvals do
    for r2 in Rvals do
        u := Lump(k2, l2, r2, 0):
        P := plot3d(u(x, y), x = -15 .. 15, y = -15 .. 15, color = gold,
                    title = cat("Lump: k[2]=", k2, ", l[2]=", l2, ", r[2]=", r2),
                    labels = ["x", "y", "u(x,y)"], axes = boxed, grid = [60, 60]):
        plotsList := [op(plotsList), P]:
    od:
  od:
od:

display(plotsList, insequence = false);

# Stap 1: Stel parameterwaarden in
k2_val := 1: l2_val := 1: r2_val := 0.5: tval := 0:

# Stap 2: Genereer lump-oplossing
ulump := Lump(k2_val, l2_val, r2_val, tval):

# Stap 3: Originele PDE
pde := diff(diff(u(x, y, z, t), t) + 6*u(x, y, z, t)*diff(u(x, y, z, t), x) + diff(u(x, y, z, t), x $ 3), x)
       - lambda*diff(u(x, y, z, t), y $ 2) + diff(alpha*diff(u(x, y, z, t), x) + beta*diff(u(x, y, z, t), y) + delta*diff(u(x, y, z, t), z), x):

# Stap 4: Substitueer lump en evalueer op z = 0
subcheck := eval(pde, u(x, y, z, t) = ulump(x, y)):
subcheck := eval(subcheck, z = 0):

# Stap 5: Controleer of het resultaat ~0
simplify(subcheck);

restart:
with(PDEtools): with(plots):

# ✅ 1. Lump-oplossing als functie
Lump := proc(k2, l2, r2, tval)
    local G, fx, fxx, u;
    G := exp(-(1/k2)*(tval*k2^4 + tval*k2^2 + tval*k2*l2 + tval*k2*r2 - tval*l2^2 - x*k2^2 - y*k2*l2)) +
         exp(-1) +
         exp(-(1/k2)*(tval*k2^4 + tval*k2^2 + tval*k2*l2 + tval*k2*r2 - tval*l2^2 - x*k2^2 - y*k2*l2)):
    fx := diff(G, x): fxx := diff(G, x$2):
    u := unapply(2*fxx/G - 2*(fx/G)^2, x, y):
    return u:
end proc:

# �� 2. Originele PDE
pde := diff(diff(u(x, y, z, t), t) + 6*u(x, y, z, t)*diff(u(x, y, z, t), x) + diff(u(x, y, z, t), x $ 3), x)
       - lambda*diff(u(x, y, z, t), y $ 2) + diff(alpha*diff(u(x, y, z, t), x) + beta*diff(u(x, y, z, t), y) + delta*diff(u(x, y, z, t), z), x):

# �� 3. Vaste lump parameters
k2_val := 1: l2_val := 1: r2_val := 0.5: tval := 0:
ulump := Lump(k2_val, l2_val, r2_val, tval):

# �� 4. Testen over alpha, beta, lambda
Avals := [0, 1, 2]:
Bvals := [0, 1, 2]:
Lvals := [0, 1, 2]:
plotsList := []:

for a in Avals do
  for b in Bvals do
    for lam in Lvals do
        alpha := a: beta := b: lambda := lam:
        check := eval(pde, u(x, y, z, t) = ulump(x, y)):
        check := eval(check, z = 0):
        check := simplify(check):
        P := plot3d(check, x = -5 .. 5, y = -5 .. 5,
                    title = cat("Residue voor α=", a, ", β=", b, ", λ=", lam),
                    axes = boxed, labels = ["x", "y", "Residue"]):
        plotsList := [op(plotsList), P]:
    od:
  od:
od:

# �� 5. Alles tonen
display(plotsList, insequence = false, title = "Residuanalyse van Lump-substitutie");

restart;
with(plots):

# Parameters voor de lumps
k1 := 1: l1 := 1: r1 := 0.5: eta1 := 0:
k2 := 1.2: l2 := -1: r2 := 0.3: eta2 := 0:
t := 0:

# Definieer de fasen (theta1 en theta2)
theta1 := t*k1^4 + t*k1^2 + t*k1*l1 + t*k1*r1 - t*l1^2 + x*k1^2 + y*k1*l1 + eta1:
theta2 := t*k2^4 + t*k2^2 + t*k2*l2 + t*k2*r2 - t*l2^2 + x*k2^2 + y*k2*l2 + eta2:

# F-functie van Hirota (twee lumps + interactie)
b := 1: # Interactiecoëfficiënt
F := 1 + exp(theta1) + exp(theta2) + b*exp(theta1 + theta2):

# Lump-oplossing via Hirota's formule
fx := diff(F, x):
fxx := diff(F, x$2):
u := unapply(2*fxx/F - 2*(fx/F)^2, x, y):

# 3D-plot
P3D := plot3d(u(x, y), x = -10 .. 10, y = -10 .. 10,
              title = "3D-plot van tweelump-oplossing",
              axes = boxed, grid = [60, 60], color = gold,
              labels = ["x", "y", "u(x,y)"]);

# Contourplot
Pcontour := contourplot(u(x, y), x = -10 .. 10, y = -10 .. 10,
                        contours = 20, coloring = [blue, white, red],
                        title = "Contourplot van tweelump-oplossing",
                        axes = boxed);

# Toon beide plots
#display([P3D, Pcontour], title = "2-Lump Solution Visualisatie");

Error, (in plots:-display) cannot display 2-D and 3-D plots together

restart:
with(PDEtools):

# Zorg dat f als functie wordt herkend
declare(f(x,y,z,t), 'function'):

# �� Hoofdroutine: controleer of een f(x,y,z,t) voldoet aan Hirota-bilineaire vorm
CheckHirotaForm := proc(f_expr::algebraic, vars::list)
    description "Onderzoekt of een functie f(x,y,z,t) voldoet aan de Hirota-bilineaire PDE";

    local x, y, z, t, f;
    local A, B, A1, B1, C1;
    local D_op, Dff, BilinearPDE, Result;

    # Variabelen extraheren
    (x, y, z, t) := op(vars):

    # Zet expressie om naar functionele vorm
    f := unapply(f_expr, x, y, z, t):

    # Stel Hirota-operator op voor f*f
    D_op := A*convert(Diff(f(x,y,z,t)*f(x,y,z,t), x,x,x,x), D) +
            A1*convert(Diff(f(x,y,z,t)*f(x,y,z,t), x, y), D) +
            B1*convert(Diff(f(x,y,z,t)*f(x,y,z,t), x, z), D) +
            C1*convert(Diff(f(x,y,z,t)*f(x,y,z,t), t, x), D) +
            convert(Diff(f(x,y,z,t)*f(x,y,z,t), t, t), D) +
            convert(Diff(f(x,y,z,t)*f(x,y,z,t), x, x), D):

    # Stel volledige Hirota PDE op
    BilinearPDE := convert(Diff(3*A*D_op / (B*f(x,y,z,t)^2), x$2), D):

    # Vereenvoudig resultaat
    Result := simplify(BilinearPDE):

    # Toon uitkomst
    print("Bilineaire vergelijking na substitutie:");
    print(Result);

    # Controleer of nul
    if Result = 0 then
        print("✅ Geldige oplossingsfunctie!");
    else
        print("❌ Niet een exacte oplossing.");
    end if;

    return Result;
end proc:

# ✨ Subprocedure: test klassieke 1-soliton oplossing
Test1Soliton := proc()
    description "Genereert en test een klassieke 1-soliton-ansatz";

    local f_expr;

    # Stel 1-soliton-ansatz op
    f_expr := 1 + exp(k*x + l*y + m*z + omega*t + delta):

    print("➡️ Test van 1-soliton functie: f(x,y,z,t) = 1 + exp(kx + ly + mz + ωt + δ):");
    return CheckHirotaForm(f_expr, [x, y, z, t]);
end proc:

# Gebruik:
# Test1Soliton();

# Stel waarden in voor de parameters (optioneel symbolisch of numeriek)
A := 1: B := 6:
A1 := 0: B1 := 0: C1 := 0:

# Voer de test uit
Test1Soliton();