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

# Parameters
assume(a::real, b::real, alpha::real, beta::real, delta::real, mu::real);
assume(k[1]::real, k[2]::real, r[1]::real, r[2]::real, l[1]::real, l[2]::real);
assume(x::real, y::real, z::real, t::real);

a := 1:   b := 0.5:
alpha := 1:  beta := 0.8: delta := 0.6: mu := 1:
k[1] := 1:  k[2] := 2:
r[1] := 0.3:  r[2] := 0.7:

# Correct l[1] and l[2] to ensure real solutions
l[1] := -3:
l[2] := -7:

# Soliton solution components
thetai := k[i]*(t*w[i] + y*l[i] + z*r[i] + x);

# Compute w[i]
eqw := w[i] = (-1 + sqrt(-4*beta*mu*l[i] -4*delta*mu*r[i] -4*mu*k[i]^2 -4*alpha*mu +1))/(2*mu);
w[1] := eval(rhs(eqw), i=1);
w[2] := eval(rhs(eqw), i=2);

# Interaction term Bij
Bij := (i, j) -> simplify(-24*mu/(sqrt(1+(-4*beta*l[j]-4*delta*r[j]-4*alpha)*mu)*sqrt(1+(-4*beta*l[i]-4*delta*r[i]-4*alpha)*mu) -1 + ((2*r[i]+2*r[j])*delta + (2*l[i]+2*l[j])*beta +4*alpha)*mu), size);

# Theta functions
theta1 := eval(thetai, i=1);
theta2 := eval(thetai, i=2);

# Define f
eqf := theta1*theta2 + Bij(1, 2);

# Compute u
f := eqf;
f_x := diff(f, x);
f_xx := diff(f_x, x);
u_expr := (2*f_xx/f - 2*(f_x)^2/f^2);

# Simplify
equ := simplify(convert(u_expr, rational), size);

# Plot 3D solution
p1 := plot3d(eval(equ, [z=0, y=0]), x=-10..10, t=-10..10,
       axes=boxed, grid=[60,60], shading=zhue,
       
       labels=["x", "t", "u(x,t)"], title="Soliton Solution in 3D");

# 2D plot at t=0
p2 := plot(eval(equ, [z=0, y=0, t=0]), x=-10..10,
       color=blue, labels=["x", "u(x,0)"],
       title="Soliton Profile at t=0");

# Contour plot
p3 := contourplot(eval(equ, [z=0, y=0]), x=-10..10, t=-10..10,
       grid=[60,60], coloring=[blue, red],
       labels=["x", "t"],
       title="Contour Plot of Soliton");
 

The goal is to implement an expression of a discrete scalar plane-wave solution to the Klein-Gordon equation that makes use of the creation and annihilation operators in the Physics package.

We don't bother solving the Klein-Gordon equation explicitly. Instead, we just take a solution from literature. This approach is sufficient to show the difficulty in defining the the creation and annihilation operators.

restart;
with(PDEtools):
with(plots):
declare(u(x, y, t), quiet);

# Parameters
a[1] := -1:
b[1] := -2:
a[2] := -2/3:
b[2] := -1:
alpha := 1:
beta := 1:

# Definitie van de 2-lump oplossing
u_expr := -2*((9/65 - 18*I/65)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + (9/65 + 18*I/65)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + (7933424/1891125 + 36508832*I/1891125)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x) + (6/65 + 9*I/65)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + (6/65 - 9*I/65)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x) + (8059833/1092650 - 84073113*I/2185300)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + (7933424/1891125 - 36508832*I/1891125)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x) + (8059833/1092650 + 84073113*I/2185300)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x))/(9*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x)/65 + (27783/1625 - 11619*I/1625)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + 135*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x)/52 + (74331/109265 + 118467*I/109265)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x) + (74331/109265 - 118467*I/109265)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*((-2/3 + I)*t + (5/9 + 4*I/3)*y + t + (2/3 - I)*x) + 26*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*((-1 + 2*I)*t + (3 + 4*I)*y + t + (1 - 2*I)*x)/45 + (27783/1625 + 11619*I/1625)*(-(1 + 2*I)*t + (3 - 4*I)*y + t + (1 + 2*I)*x)*(-(2/3 + I)*t + (5/9 - 4*I/3)*y + t + (2/3 + I)*x) + 126229597/50430):

# Trajectorieën
eq1 := y = -beta/((a[1]^2 + b[1]^2)*alpha + 2*a[1]*beta)*(x - (sqrt(3)*sqrt(-a[1]*beta*(a[1]^2 + b[1]^2)^3)/(beta*(a[1]^2 + b[1]^2)*b[1]))):
eq2 := y = -beta/((a[2]^2 + b[2]^2)*alpha + 2*a[2]*beta)*(x - (sqrt(3)*sqrt(-a[2]*beta*(a[2]^2 + b[2]^2)^3)/(beta*(a[2]^2 + b[2]^2)*b[2]))):

MultiLumpPlot2 := proc(tijden::list, kleuren::list, x_range::list, y_range::list, custom_grid::list := [100, 100])
    local contourplots, tijd, i, trajectory_plot1, trajectory_plot2, all_plots, x_min, x_max, y_min, y_max;
    
    x_min, x_max := op(x_range);
    y_min, y_max := op(y_range);
    
    contourplots := [seq(contourplot(eval(u_expr, t = -tijden[i]), x = x_min .. x_max, y = y_min .. y_max, contours = 30, color = kleuren[i], grid = custom_grid, transparency = 0.1), i = 1 .. nops(tijden))];

    trajectory_plot1 := implicitplot(eq1, x = x_min .. x_max, y = y_min .. y_max, color = magenta, thickness = 2);
    trajectory_plot2 := implicitplot(eq2, x = x_min .. x_max, y = y_min .. y_max, color = black, thickness = 2);
    
    return display(op(contourplots), trajectory_plot1, trajectory_plot2, labels = ["x", "y"], scaling = constrained, size = [1200, 800]);
end proc:

MultiLumpPlot2([300], ["blue"], [-200, 200], [-100, 100]);

MultiLumpPlot2([100, 200, 300,500, -100, -200,-300 ], ["red", "green", "blue", "black", "red", "green", "blue"], [-350, 350], [-350, 350]);#