janhardo

a procedure for F_2=theta[1]*...

Answered: janhardo 545

April 20 2025

0 0

this is probably where you are looking for ?

GetF_N := proc(N)

>

(1)

>

(2)

>

(3)

>

(4)

>

(5)

>

(6)

Change b[1,2]*b[2,3]*b[1,3] to b[1,2,3], etc.

>	evalindets(f,``, proc(term) local inds:=indets(term,specindex(A)); A[op(op~(inds))]eval(term,inds=~1); end proc);

(7)

>

restart:

# Number of components
N := 2:

# Define ε-dependent variables
epsilon := 'epsilon':
for i from 1 to N do
beta[i] := epsilon*theta[i];
end do:
B[1,2] := epsilon^2*Lambda12:

# Define the general exponent
f0 := exp(add(mu[i]*beta[i], i = 1 .. N) +
add(add(mu[i]*mu[j]*B[i,j], i = 1 .. j-1), j = 1 .. N)):

# Generate all binary combinations of mu values
mus := Iterator:-BinaryGrayCode(N):
muvals := seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus):

# Sum over all binary cases
f := add(eval(f0, muval), muval in muvals):

# Expand as Taylor series in ε and extract ε² term (F2)
f_series := series(f, epsilon, 3):
F2 := simplify(coeff(f_series, epsilon, 2)); # This is F₂

(8)

>

restart:

GetF_N := proc(N)
  local i, j, epsilon, beta, theta, B, f0, f, f_series, mu, mus, muvals, F2, F2_clean, mui, muval;

  epsilon := 'epsilon':

  # Step 1: Define beta[i] = ε θ[i]
  for i from 1 to N do
    beta[i] := epsilon*theta[i]:
  end do:

  # Step 2: Define B[i,j] = ε² Λ[i,j]
  for i from 1 to N-1 do
    for j from i+1 to N do
      B[i,j] := epsilon^2*Lambda[i,j]:
    end do:
  end do:

  # Step 3: Build general Hirota term
  f0 := exp(add(mu[i]*beta[i], i = 1 .. N) +
            add(add(mu[i]*mu[j]*B[i,j], i = 1 .. j-1), j = 1 .. N)):

  # Step 4: Enumerate all binary mu combinations
  mus := Iterator:-BinaryGrayCode(N):
  muvals := seq({seq(mu[i] = mui[i], i = 1 .. N)}, mui in mus):

  # Step 5: Sum all evaluated exponential terms
  f := add(eval(f0, muval), muval in muvals):

  # Step 6: Expand in ε and extract ε² coefficient
  f_series := series(f, epsilon, 3):
  F2 := simplify(coeff(f_series, epsilon, 2)):

  # Step 7: Remove θ[i]^2 terms to get cross interactions only
  F2_clean := F2:
  for i from 1 to N do
    F2_clean := simplify(F2_clean - coeff(F2_clean, theta[i]^2)*theta[i]^2):
  end do:

  return F2_clean;
end proc:

>	GetF_N(2); # F

(9)

>

>	GetF_N(1);

(10)

>	GetF_N(3);

(11)

>

Download voorbeeld_5_N-S_worksheet_aangepast_met_procedure_20-4-2025.mw

@salim-barzani i thought i give you...

Answered: janhardo 545

April 20 2025

0 0

@salim-barzani
i thought i give you some maple code to figure this further out by yourself?

examples cpu intensive i think , hanging...

Answered: janhardo 545

April 20 2025

0 0

examples cpu intensive i think , hanging ?

>

>	# version 1

>

restart:

# Load packages

with(plots):
local gamma;
# === Parameters === #
K1 := 'K1': P1 := 'P1':
K2 := 'K2': P2 := 'P2':
K3 := 'K3': P3 := 'P3':
beta := 'beta': gamma := 'gamma': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# === Theta functions (from equation (14)) === #
theta1 := (-K1^2*gamma + P1^2*beta)/K1 + K1*y + P1*x:
theta2 := (-K2^2*gamma + P2^2*beta)/K2 + K2*y + P2*x:

# === Soliton phase η³ === #
eta3 := K3*y + P3*x:

# === B_ij terms from equation (15) === #
B13 := beta * ((K1*K3*P1*(K1 + K3)) / (K1^2*P1^2 - P3^2) + K1*K3*P1*(K1 + K3)):
B23 := beta * ((K2*K3*P2*(K2 + K3)) / (K2^2*P2^2 - P3^2) + K2*K3*P2*(K2 + K3)):

# === Construct f (Equation 13) === #
f := theta1 + theta2 + exp(eta3)*(B13*B23 + theta2*B13 + theta1*B23 + theta1*theta2):

# === Final nonlinear wave field u === #
u := simplify(6/alpha * diff(ln(f), x$2)):

# === Optional: Plot u(x,y) for given parameters === #
# You can assign values and plot:
# K1 := 1: P1 := 2: etc...
# plot3d(eval(u, [K1=1, P1=2, ...]), x = -10 .. 10, y = -10 .. 10);

(1)

>	# version 2

>

restart:

# Load necessary packages
with(plots): with(VectorCalculus):
local gamma;
# Declare symbolic parameters
K := [K1, K2, K3]: P := [P1, P2, P3]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# Define θ_i (Eq. 14): theta_i = (K_i^2 * γ + P_i^2 * β)/K_i + K_i*y + P_i*x
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:

# Define η^3 = K3*y + P3*x + ω3*t + ξ3 (simplify to symbolic form)
eta3 := K[3]*y + P[3]*x + 'omega3'*t + 'xi3':

# Define B_ij for i<j<3 (Eq. 15)
B := Matrix(3, 3):
for i from 1 to 2 do
    for j from i+1 to 3 do
        B[i,j] := beta * (
            (6*K[i]*P[i]*K[j]*((K[i]*P[j] + K[j]*P[i])))/(
            -3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# Shortcuts for readability
theta1 := theta[1]: theta2 := theta[2]:
B12 := B[1,2]: B13 := B[1,3]: B23 := B[2,3]:

# Full expression of f from Eq. (13)
f := theta1 + theta2 + exp(eta3)*(B13*B23 + theta2*B13 + theta1*B23 + theta1*theta2 + B12):

# Compute u = (6/alpha) * ∂²/∂x² ln(f)
u := simplify(6/alpha * diff(ln(f), x$2)):

(2)

Maple Code for Equation (27): 1-Lump + 2-Soliton

>

restart:

# Load packages
with(plots): with(VectorCalculus):

local gamma;
# Declare symbolic parameters
K := [K1, K2, K3, K4]: P := [P1, P2, P3, P4]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# Define θ_i terms (Eq. 14)
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:
theta1 := theta[1]: theta2 := theta[2]:

# Define η^3 and η^4
eta3 := K[3]*y + P[3]*x + 'omega3'*t + 'xi3':
eta4 := K[4]*y + P[4]*x + 'omega4'*t + 'xi4':

# Define B_ij matrix using general formula from Eq. (15)
B := Matrix(4, 4):
for i from 1 to 4 do
    for j from i+1 to 4 do
        B[i,j] := beta * (
            (6*K[i]*P[i]*K[j]*((K[i]*P[j] + K[j]*P[i])))/(
            -3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# Shortcuts for B_ij terms
B12 := B[1,2]: B13 := B[1,3]: B23 := B[2,3]:
B14 := B[1,4]: B24 := B[2,4]:

# Assume symbolic A34 (interaction coefficient)
A34 := 'A34':

# Construct f as in Equation (27)
f := theta1 + theta2 + B12
     + A34 * exp(eta3 + eta4) * (B13*B23 + B13*theta2 + B23*theta1 + theta1*theta2 + B12)
     + B13*theta2 + B12*B23 + B14*B24 + B14*theta2 + B24*theta1 + B24*theta1*theta2
     + (B13*B23 + B12) * exp(eta3)
     + (B14*B24 + B12) * exp(eta4):

# Compute u = (6/alpha) * ∂²/∂x² ln(f)
u := simplify(6/alpha * diff(ln(f), x$2)):

(3)

�� Maple Code for Equation (34): Lump + 3-Soliton Superposition

>

restart:

with(plots): with(VectorCalculus):

# === Parameters === #
K := [K1, K2, K3, K4, K5]: P := [P1, P2, P3, P4, P5]:
gamma := 'gamma': beta := 'beta': alpha := 'alpha':
x := 'x': y := 'y': t := 't':

# === θ_i terms for lump wave === #
theta := [seq((K[i]^2*gamma + P[i]^2*beta)/K[i] + K[i]*y + P[i]*x, i = 1..2)]:
theta1 := theta[1]: theta2 := theta[2]:

# === η^j soliton phases === #
eta := [seq(K[i]*y + P[i]*x + 'omega'||i*t + 'xi'||i, i=3..5)]:

# === Define B_ij terms === #

B := Matrix(5, 5):
for i from 1 to 5 do
    for j from i+1 to 5 do
        B[i, j] := beta * (
            (6*K[i]*P[i]*K[j]*(K[i]*P[j] + K[j]*P[i])) /
            (-3*P[i]^2*K[i]^2 - 2*P[i]*K[i]*beta + beta*P[i]^2 + gamma*K[i]^2)
        );
    end do;
end do:

# === A_ij interaction terms — symbolic form === #
A34 := 'A34': A35 := 'A35': A45 := 'A45':

# === Reuse B_ij notation === #
B12 := B[1,2]: B13 := B[1,3]: B14 := B[1,4]: B15 := B[1,5]:
B23 := B[2,3]: B24 := B[2,4]: B25 := B[2,5]:
B34 := B[3,4]: B35 := B[3,5]: B45 := B[4,5]:

# === Build f (Eq. 34) === #
f := theta1 + theta2 + B12
  + (B13 + theta1)*(B23 + theta2)*exp(eta[1])
  + (B14 + theta1)*(B24 + theta2)*exp(eta[2])
  + (B15 + theta1)*(B25 + theta2)*exp(eta[3])
  + A34*(B13 + theta1)*(B23 + theta2)*exp(eta[1] + eta[2])
  + A35*(B13 + theta1)*(B23 + theta2)*exp(eta[1] + eta[3])
  + A45*(B14 + theta1)*(B24 + theta2)*exp(eta[2] + eta[3])
  + A34*A35*A45*(B13 + B14 + B15)*(B23 + B24 + B25)*exp(eta[1] + eta[2] + eta[3]):

# === Final nonlinear wave field u === #
u := simplify(6/alpha * diff(ln(f), x$2)):

Download deel4_lump_waves_20-4-2025.mw

restart: # Load packages with(plot...

Answered: janhardo 545

April 20 2025

0 0

restart:

# Load packages
with(plots): with(VectorCalculus):

# Declare variables
epsilon := 'epsilon':

# Define symbolic long-wave parameters
K := [K1, K2, K3, K4]:
L := [L1, L2, L3, L4]:
R := [R1, R2, R3, R4]:
Omega := [O1, O2, O3, O4]:
X := [X1, X2, X3, X4]:

# Construct eta_i = ε Φ_i
eta := [seq(epsilon*(K[i]*x + L[i]*y + R[i]*z + Omega[i]*t + X[i]), i = 1 .. 4)]:

# A_ij matrix (interaction terms)
A := Matrix(4, 4):
for i from 1 to 4 do
    for j from i+1 to 4 do
        A[i, j] := ln(((epsilon*K[i] - epsilon*K[j])^2)/((epsilon*K[i] + epsilon*K[j])^2));
    end do;
end do:

# Power set for binary mu-values
mu := combinat:-powerset({1,2,3,4}):
f := 0:

# Construct full f-series
for S in mu do
    term := 0:
    for i in S do
        term := term + eta[i];
    end do:
    for i in S do
        for j in S do
            if i < j then
                term := term + A[i, j];
            end if;
        end do;
    end do:
    f := f + series(exp(term), epsilon, 5);  # Expand each term individually
end do:

# Simplify the whole tau-function
f_expanded := collect(simplify(f), epsilon):
f_lump := simplify(coeff(f_expanded, epsilon, 4));  # Take ε⁴ coefficient

note : code is not working yet: keeps hanging

@salim-barzani J # Param...

Answered: janhardo 545

April 20 2025

0 0

@salim-barzani
J

# Parameters
epsilon := 'epsilon':
K1 := 'K1': K2 := 'K2':
L1 := 'L1': L2 := 'L2':
R1 := 'R1': R2 := 'R2':
O1 := 'Omega1': O2 := 'Omega2':
X1 := 'X1': X2 := 'X2':

# Eta terms
Phi1 := K1*x + L1*y + R1*z + O1*t + X1:
Phi2 := K2*x + L2*y + R2*z + O2*t + X2:
eta1 := epsilon * Phi1:
eta2 := epsilon * Phi2:

# A12 term
A12 := ln(((epsilon*K1 - epsilon*K2)^2)/((epsilon*K1 + epsilon*K2)^2)):

# Exponentials (up to ε²)
exp1 := series(exp(eta1), epsilon, 3):
exp2 := series(exp(eta2), epsilon, 3):
exp12 := series(exp(eta1 + eta2 + A12), epsilon, 3):

# Total tau-function
f := 1 + exp1 + exp2 + exp12:
f_simplified := collect(simplify(f), epsilon):

# Extract ε² coefficient — this is the lump form!
f_lump := simplify(coeff(f_simplified, epsilon, 2));

@salim-barzani ...

Answered: janhardo 545

April 20 2025

0 0

@salim-barzani

If you think this is not what , then it ...

Answered: janhardo 545

April 20 2025

0 0

If you think this is not what i am looking for , then it can be deleted of course

...

Answered: janhardo 545

April 20 2025

0 0

>

restart:

with(plots): with(VectorCalculus):

# Define complex wave parameters for two lump pairs
K := [1 + I, 1 - I, 2 + 2*I, 2 - 2*I]:
L := [-1/2 + 2*I, -1/2 - 2*I, 1 + I, 1 - I]:
R := [-2 + I, -2 - I, 0 + I, 0 - I]:
Omega := [1, 1, 1, 1]:
local Chi ;
Chi := [0, 0, 0, 0]: # No initial phase

# Define eta_i terms using ε substitution
epsilon := 'epsilon':
eta := [seq(epsilon*(K[i]*x + L[i]*y + R[i]*z + Omega[i]*t + Chi[i]), i=1..4)]:

# Define A_ij terms (interaction phase shifts)
A := Matrix(4, 4):
for i from 1 to 4 do
    for j from 1 to 4 do
        if i < j then
            ki := epsilon*K[i]: kj := epsilon*K[j]:
            A[i, j] := ln(((ki - kj)^2)/((ki + kj)^2));
        end if;
    end do;
end do:

# Define the tau-function f with binary mu-summation over 4 variables
mu := combinat:-powerset({1,2,3,4}):
f := 0:
for S in mu do
    term := 0:
    for i in S do
        term := term + eta[i]:
    end do:
    for i in S do
        for j in S do
            if i < j then
                term := term + A[i, j];
            end if;
        end do;
    end do:
    f := f + exp(term);
end do:

# Expand in ε up to 5th order
f_exp := series(f, epsilon, 5):
f_lump := simplify(coeff(f_exp, epsilon, 4)): # Take ε⁴ coefficient

# u(x,y) = 2nd x-derivative of ln(f)
u := simplify(2 * diff(ln(f_lump), x$2)):

# Evaluate at z=0, t=0 to visualize
u_plot := subs(z=0, t=0, u):

# 3D plot
plot3d(u_plot, x = -6 .. 6, y = -6 .. 6,
       axes = boxed, labels = ["x", "y", "u(x,y)"],
       title = "Two-Lump Solution");

>

restart:

with(plots): with(VectorCalculus):

# Define complex wave parameters for two lump pairs
K := [1 + I, 1 - I, 2 + 2*I, 2 - 2*I]:
L := [-1/2 + 2*I, -1/2 - 2*I, 1 + I, 1 - I]:
R := [-2 + I, -2 - I, 0 + I, 0 - I]:
Omega := [1, 1, 1, 1]:
Chi_ := [0, 0, 0, 0]: # Renamed from Chi to Chi_

# Define eta_i terms using ε substitution
epsilon := 'epsilon':
eta := [seq(epsilon*(K[i]*x + L[i]*y + R[i]*z + Omega[i]*t + Chi_[i]), i=1..4)]:

# Define A_ij terms (interaction phase shifts)
A := Matrix(4, 4):
for i from 1 to 4 do
    for j from 1 to 4 do
        if i < j then
            ki := epsilon*K[i]: kj := epsilon*K[j]:
            A[i, j] := ln(((ki - kj)^2)/((ki + kj)^2));
        end if;
    end do;
end do:

# Define the tau-function f with binary mu-summation over 4 variables
mu := combinat:-powerset({1,2,3,4}):
f := 0:
for S in mu do
    term := 0:
    for i in S do
        term := term + eta[i]:
    end do:
    for i in S do
        for j in S do
            if i < j then
                term := term + A[i, j];
            end if;
        end do;
    end do:
    f := f + exp(term);
end do:

# Expand in ε up to 5th order
f_exp := series(f, epsilon, 5):
f_lump := simplify(coeff(f_exp, epsilon, 4)): # Take ε⁴ coefficient

# u(x,y) = 2nd x-derivative of ln(f)
u := simplify(2 * diff(ln(f_lump), x$2)):

# Evaluate at z=0, t=0 to visualize
u_plot := subs(z=0, t=0, u):

# 3D plot
plot3d(u_plot, x = -6 .. 6, y = -6 .. 6,
       axes = boxed, labels = ["x", "y", "u(x,y)"],
       title = "Two-Lump Solution");

>

Download deel3_vraag_beantwoord_Two_lump_solution_20-4-2025.mw

@salim-barzani Is this better ..?&n...

Answered: janhardo 545

April 20 2025

0 0

@salim-barzani
Is this better ..?

>

restart;

# Parameters
epsilon := 'epsilon':
K1 := 'K1': L1 := 'L1': R1 := 'R1': Omega1 := 'Omega1': X1 := 'X1':
K2 := 'K2': L2 := 'L2': R2 := 'R2': Omega2 := 'Omega2': X2 := 'X2':

# Substitutions
k1 := epsilon*K1: l1 := epsilon*L1: r1 := epsilon*R1: omega1 := epsilon*Omega1: chi1 := epsilon*X1:
k2 := epsilon*K2: l2 := epsilon*L2: r2 := epsilon*R2: omega2 := epsilon*Omega2: chi2 := epsilon*X2:

# Eta expressions
eta1 := k1*x + l1*y + r1*z + omega1*t + chi1:
eta2 := k2*x + l2*y + r2*z + omega2*t + chi2:

# A12 term (second-order small term)
A12 := ln(((k1 - k2)^2)/(k1 + k2)^2):

# Full exponential sum (Eq. 6 reduced to N=2)
f := 1 + exp(eta1) + exp(eta2) + exp(eta1 + eta2 + A12):

# Series expansion as epsilon -> 0
f_exp := series(f, epsilon, 3):

# Extract order epsilon^2 term (leading order for lump wave)
simplify(f_exp, symbolic);

series(3+(K1-K2)^2/(K1+K2)^2+(2*(K1^2+K2^2)*(K1*x+K2*x+(Omega1+Omega2)*t+(L1+L2)*y+(R2+R1)*z+X1+X2)/(K1+K2)^2)*epsilon+((1/2)*(K1*x+L1*y+Omega1*t+R1*z+X1)^2+(1/2)*(K2*x+L2*y+Omega2*t+R2*z+X2)^2+(1/2)*((Omega1+Omega2)*t+(K1+K2)*x+(L1+L2)*y+(R2+R1)*z+X1+X2)^2*(K1-K2)^2/(K1+K2)^2)*epsilon^2+O(epsilon^3),epsilon,3)

(1)

>

restart:

# Complexe parameters uit Fig. 4
K1 := 1 + I: L1 := -1/2 + 2*I: R1 := -2 + I:
K2 := conjugate(K1): L2 := conjugate(L1): R2 := conjugate(R1):

# Verplaatsingssnelheden, kies symbolisch (voor illustratie)
Omega1 := 'Omega1': Omega2 := 'Omega2':
X1 := 0: X2 := 0: # chi = 0 in Fig. 4

# Maak phi1 en phi2 expliciet
phi1 := K1*x + L1*y + R1*z + Omega1*t + X1:
phi2 := K2*x + L2*y + R2*z + Omega2*t + X2:

# Bereken product
phi_product := expand(phi1 * phi2):

# Bepaal A12 volgens standaard Hirota-formule (voor KP-type systeem)
# Merk op: met k_i = epsilon*K_i → we gebruiken limietvorm van A_ij
k1 := epsilon*K1: k2 := epsilon*K2:
A12 := ln( ((k1 - k2)^2) / (k1 + k2)^2 ):
A12_series := series(A12, epsilon, 3): # expanderen tot ε²
Lambda12 := simplify(limit(A12_series, epsilon = 0)): # ε²-termen bewaren als constante

# Print de resultaten
phi1, phi2, phi_product, Lambda12;

(1+I)*x+(-1/2+2*I)*y+Omega1*t+(-2+I)*z, (1-I)*x+(-1/2-2*I)*y+Omega2*t+(-2-I)*z, -2*z*Omega2*t+Omega1*t*x+x*Omega2*t-(1/2)*y*Omega2*t+Omega1*t^2*Omega2-(1/2)*Omega1*t*y-2*Omega1*t*z+2*x^2+3*x*y-2*x*z+(17/4)*y^2+6*y*z+5*z^2-I*Omega1*t*z+(2*I)*y*Omega2*t+I*z*Omega2*t+I*x*Omega2*t-I*Omega1*t*x-(2*I)*Omega1*t*y, I*Pi

(2)

>

restart:

with(plots): with(VectorCalculus):

# Complexe parameters
K1 := 1 + I: L1 := -1/2 + 2*I: R1 := -2 + I:
K2 := conjugate(K1): L2 := conjugate(L1): R2 := conjugate(R1):

# Symbolische Omega waarden
Omega1 := 1: Omega2 := 1:

# Chi = 0
X1 := 0: X2 := 0:

# φ1 en φ2
phi1 := K1*x + L1*y + R1*z + Omega1*t + X1:
phi2 := K2*x + L2*y + R2*z + Omega2*t + X2:

# f(x, y, z, t)
f := simplify(phi1 * phi2 + 1): # Lambda12 = 1 aangenomen voor visualisatie

# u(x, y, z, t) = 2nd x-derivative of ln(f)
u := simplify(2 * diff(ln(f), x$2)):

# Fixeer z en t zodat u(x, y) plotbaar is
u_plot := subs(z=0, t=0, u):

# Plot op domein
plot3d(u_plot, x = -5 .. 5, y = -5 .. 5,
axes = boxed, labels = ["x", "y", "u(x,y)"],
title = "One-Lump Solution");

>

Download deel_2_vraag_beantwoorrden_20-4-2025.mw

@salim-barzani it's ai style , ...

Answered: janhardo 545

April 19 2025

0 0

@salim-barzani
it's ai style , so ?
i tried different versions to get for r2 , as close as the paper
Probably it can be better , but how you can derive r2 is clear

>

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.

r1+(1/4)*3^(1/2)*(-(k1-k2)^2*a^4+6*b^2*(k1-k2)^2*a^2-(k1-k2)^2*b^4+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2)

r1-(1/4)*3^(1/2)*(-(k1-k2)^2*a^4+6*b^2*(k1-k2)^2*a^2-(k1-k2)^2*b^4+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2)

[r1+(1/4)*3^(1/2)*(-(k1-k2)^2*a^4+6*b^2*(k1-k2)^2*a^2-(k1-k2)^2*b^4+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2), r1-(1/4)*3^(1/2)*(-(k1-k2)^2*a^4+6*b^2*(k1-k2)^2*a^2-(k1-k2)^2*b^4+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2)]

[r1+(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k1-k2)^2+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2), r1-(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k1-k2)^2+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2)]

(1)

>

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) ];

[(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k2-k1)^2*lambda+4*lambda*(r1+s1-s2))^(1/2)/lambda, -(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k2-k1)^2*lambda+4*lambda*(r1+s1-s2))^(1/2)/lambda]

(2)

>

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);

r1+(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k2-k1)^2+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2), r1-(1/4)*3^(1/2)*(-(a^4-6*a^2*b^2+b^4)*(k2-k1)^2+4*r1+4*s1-4*s2)^(1/2)/lambda^(1/2)

(3)

Download berekinhg_mprimes_vraag_formule_bereekn_in_pde_engelse_versie19-14-2025.mw

@salim-barzani ...

Answered: janhardo 545

April 19 2025

0 0

@salim-barzani

a equation for A12 = 0 , and solve this&...

Answered: janhardo 545

April 19 2025

0 0

a equation for A12 = 0 , and solve this

assume : Lambda > 0

> ....

Commented: janhardo 545

April 19 2025

>

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);

(1)

>

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");