Korteweg-de Vries Equation Soliton Solution

Abstract

Matlab code was developed to solve the Korteweg-de Vries equation for a soliton solution using an iterative method based on the two-dimensional Taylor expansion of the solution. Equations of this form frequently appear in physics, making their solutions of great importance. The code divides the solution into displacement grid points, evaluates the disturbance's position at discrete time intervals, and assesses its stability based on user-defined parameters. The results are compared to theoretical predictions, and the method can be applied to any soliton solution of this differential equation, provided the problem is appropriately scaled.

1. Introduction

The objective of this code is to solve the Korteweg-de Vries (KdV) equation (1) for a soliton solution and examine its evolution over time to determine its stability:

$$frac{partial u}{partial t} + u frac{partial u}{partial x} + delta^2 frac{partial^3 u}{partial x^3} = 0$$ (1)

Soliton solutions are wave-pulse solutions of non-linear wave equations. They satisfy the KdV equation (1), which is a reformulated version of the general equation describing a disturbance of height $$eta(r, t)$$ in a one-dimensional channel:

$$frac{partial eta}{partial t} + c frac{partial eta}{partial r} + alpha eta frac{partial eta}{partial r} + delta^2 frac{partial^3 eta}{partial r^3} = 0$$ (2)

Here, we use a change of variables, $$x = r - ct$$ and $$u = alpha eta$$, to simplify the equation.

The KdV equation describes the behavior of small but finite amplitude waves that are asymptotic over a long time period. It finds applications in various areas of physics, including shallow water waves and magnetohydrodynamic waves in collisionless plasma.

In this code, we focus on KdV solutions that are asymptotically constant at infinity, $$u_infty$$ at $$x = pm infty$$:

$$u = u_infty + (u_0 - u_infty) text{sech}^2 left(frac{x - x_0}{sigma}right)$$ (3)

For simplicity, we set $$u_infty = 0$$, and some properties of this solution are:

$$sigma = delta left(frac{u_0 - u_infty}{12}right)^{-frac{1}{2}}$$ (4)

$$c = u_infty + frac{u_0 - u_infty}{3}$$ (5)

The parameter $$sigma$$ determines the width of the soliton peak, and the speed $$c$$ is linearly proportional to its amplitude.

Additionally, the code explores phenomena like 'Classical Overtaking,' where the first two terms of the KdV equation dominate for small $$delta$$, causing regions of steeper gradients in $$u$$, and the third term prevents discontinuity formation.

The interaction of solitons, including their linear speed exchange during overlap, is not tested in this code but is important in the system's behavior.

2. Methods - Solving the Problem Numerically

To numerically solve the Korteweg-de Vries equation (1), an integration method was established with equal time and space intervals denoted as ∆t and ∆x. This discretizes the solution, evaluating it only at specific grid points: $$x = nDelta x$$ and $$t = vDelta t$$, where:

$$u_{vn} = u(nDelta x, vDelta t)$$ (6)

By Taylor expanding this two-variable function, we obtain the following approximations to substitute into the KdV equation:

$$frac{partial u}{partial t} approx u_{v+1, n} - u_{v-1, n} frac{2Delta t}{2}$$ (7)

$$frac{partial u}{partial x} approx frac{u_{vn+1} - u_{vn-1}}{2Delta x}$$ (8)

$$u approx u_{vn+1} + u_{vn} + u_{vn-1} frac{1}{3}$$ (9)

$$frac{partial^3 u}{partial x^3} approx frac{u_{vn+2} - 2u_{vn+1} + 2u_{vn-1} - u_{vn-2}}{2(Delta x)^3}$$ (10)

These approximations are not unique (see Appendix A for their derivation) but are chosen for stability over multiple iterations. Substituting these into the KdV equation yields an iterative solution that allows us to numerically find solutions:

$$u_{v+1, n} approx u_{v-1, n} - frac{2Delta t}{2Delta x}left(frac{1}{3}(u_{vn+1} + u_{vn} + u_{vn-1})(u_{vn+1} - u_{vn-1}) - delta^2frac{1}{2(Delta x)^3}(u_{vn+2} - 2u_{vn+1} + 2u_{vn-1} - u_{vn-2})right)$$ (11)

However, a problem arises with this method (equation 11) as it requires displacement points $$x$$ that are two positions below and two positions above the current point. If we only have a solution from $$x = x_0$$ to $$x = x_0 + nDelta x$$, calculating the terms at $$x = x_0$$, $$x = x_0 + Delta x$$, $$x = x_0 + (n-1)Delta x$$, and $$x = x_0 + nDelta x$$ at later times becomes problematic.

To address this issue, periodic boundary conditions are applied in this solution, where $$u_{v1} = u_{vN+1}$$, with $$N$$ grid points. From a practical perspective, this means that the vector describing the wave at a specific time must have a length of $$N + 4$$ during iterations and then be reduced to $$N$$ afterward for plotting purposes.

Another aspect to address is that the iterative method to find the next time point usually requires the current time point and the one immediately before it. However, this is not possible in the initial time instance, so a slight modification is necessary. The first substitution changes as follows (7):

$$frac{partial u}{partial t} approx u_{v+1, n} - u_{vn} Delta t$$ (12)

This modification involves altering the equation (11) by replacing $$2Delta t$$ with $$Delta t$$ and substituting $$u_{v-1, n}$$ with $$u_{vn}$$ as the first term.

These adjustments account for the unique conditions at the initial time step. The equations implemented in the code are derived from this theoretical foundation.

A dedicated function, "SolveKdV" (see Appendix B), was created to take initial wave values as inputs and perform the necessary iterations to compute the solution at later times. To expedite the computational process, a float variable is used, reducing the number of values processed. This variable contains only the current two time steps required for the iteration and the current time itself. Instead of swapping entire columns, indexing is rotated so that older values are replaced, and the two columns needed for the next iteration are utilized.

After the first iteration, the value of the first index is set to the second, the second to the third, and the third to the first. This cycle repeats for each time step, calculating the next set of values. The results are stored in a matrix where each column represents a time step.

Additionally, a script named "ScriptAnimateSolutionKdV" (see Appendix C) was developed to execute this function with user-defined parameters, including the number of grid points for solution division ($$N$$), the time interval for generating wave images ($$∆t$$), and the dispersion constant ($$δ$$). This script computes the initial solution based on the provided parameters and inputs it into the function. It then produces an animation of the wave at various time intervals, allowing users to observe its evolution over time.

The integration method used to solve the KdV equation has proven to be effective and stable when the time step is appropriately chosen for the desired calculation duration. It is important to strike a balance between solution stability and computation time. Notably, the wave's behavior aligns with expectations as the peak moves to the right with a speed of _{c ≈ 1/3} (determined by taking the ratio of displacement over time), which agrees with the theoretical prediction ⁽⁵⁾.

The wave's shape remains relatively constant throughout the simulation, and the final graph illustrates the periodic boundary conditions enforced, with the wave reappearing on the left side as it exits on the right.

Various parameters impact the stability of the solution, and they interact simultaneously. The iterative formula (11) outlines their interactions and potential sources of error, requiring verification of their expected behavior. One crucial parameter is the time step, as an excessively large time step destabilizes the solution. The code was run with _{Δt = 0.002}, demonstrating stability for over two seconds, while with _{Δt = 0.005}, the solution rapidly deteriorated. To maintain stability, a small time step is necessary, albeit at the cost of longer computation time for late times.

The number of grid points at which the wave is evaluated is another parameter that can be adjusted. Smaller _Δx values inherently introduce more error since the iterative formula (11) involves dividing by _Δx. Increasing the number of grid points from _{N = 50} to _{N = 75} led to faster instability, indicating that a balance must be struck between resolution and stability.

Having more grid points enhances solution resolution, as evident when comparing the peak visibility for _{N = 25} and _{N = 75}, but this increased resolution comes at the expense of stability. Achieving a high-resolution wave may not always be feasible, depending on the interplay of other parameters.

The impact of varying the value of δ was also explored, as it significantly affects the stability of the solution. Notably, δ is the one physical variable that is explicitly defined in the equation. Changing the number of grid points or the time step can be adjusted regardless of the situation, but δ is determined by the medium. Therefore, understanding how other parameters might need to be adjusted to accommodate different dispersion constants is useful. The width of the peak is inversely proportional to the value of δ, as predicted by theory.

If the width becomes larger than the displacement axis, the iterative formula fails to work due to the periodic boundary conditions, resulting in a discontinuity in the gradient of the solution. If δ is too large, either the axis can be extended, or the wave's amplitude can be changed to compensate.

While a lower value of the dispersion constant yields a more stable solution, the narrower peak requires higher resolution for proper visualization. This suggests that the problem should be scaled to keep δ in the range of 0.02-0.04, ensuring that parameters like _{N = 50} and _{∆t = 0.002} work optimally at this scale.

Another predicted aspect of the solution that can be tested is the phenomenon of "Classical Overtaking." This can be observed by setting δ = 0 and manually adjusting the amplitude within the code.

As predicted, in regions with a negative gradient, the gradient becomes steeper, and later, the wave begins to split into multiple waves with decreasing amplitude. The solution appears to remain stable, even though the effect is unusual. This suggests that it is the third term in equation (1) that maintains the wave in roughly its original shape. However, the same third term is also responsible for the instability of the solution when present.

Taking all these factors into account, it is evident that this method is quite effective for solving the Korteweg-de Vries equation. However, there is a trade-off between wave resolution and stability at later times, with resolution primarily defined by the dispersion constant and the smallest allowable time step. Increasing the number of grid points, which introduces instability, can be offset by reducing the time step. Yet, this may lead to longer computation times. The parameters discussed in this section can be loosely applied to any KdV equation that falls within the range of δ between 0.02 and 0.04. While the number of grid points can be adjusted slightly from _{N = 50} and the time step from _{∆t = 0.002}, these values provide a stable solution over an extended time period.

There are alternative methods for solving the KdV equation that can be more efficient but often require advanced mathematics. One such alternative is an exact solution using the Bäcklund transform, as outlined in 'The Korteweg-de Vries Equation: History, exact Solutions, and graphical Representation' by Klaus Brauer, University of Osnabrück/Germany [3]. Both this alternative method and the one presented here can be employed to study how solitons interact, although this aspect has not been explored in this study.

A. Derivation of the Partial Differential Substitutions

The general 2-D Taylor expansion up to the third order is

_{u(x + Δx, t + Δt) ≈ u(x, t) +}

^{(Δx ∂u/∂x + Δt ∂u/∂t) +}

^{1/2 (Δx² ∂²u/∂x² + 2ΔxΔt ∂²u/∂x∂t + Δt² ∂²u/∂t²) +}

^{1/6 (Δx³ ∂³u/∂x³ + 3Δx²Δt ∂³u/∂x²∂t + 3ΔxΔt² ∂³u/∂x∂t² + Δt³ ∂³u/∂t³) + ...}

(13)

This can be used to derive all of the partial substitutions, equations 7-10. First, expanding only in +Δt and -Δt

_{u(x, t + Δt) ≈ u(x, t) + Δt∂u/∂t}

(14)

_{u(x, t - Δt) ≈ u(x, t) - Δt∂u/∂t}

(15)

and subtracting (15) from (14) will yield the substitution (7). Doing exactly the same thing but with Δx instead

will yield (8).

The third partial derivative can be expressed as a combination of the following:

_{u(x + 2Δx, t) ≈ u(x, t) + 2Δx ∂u/∂x + 2(Δx)² ∂²u/∂x² + 4/3 (Δx)³ ∂³u/∂x³}

(16)

_{u(x + Δx, t) ≈ u(x, t) + Δx ∂u/∂x + 1/2 (Δx)² ∂²u/∂x² + 1/6 (Δx)³ ∂³u/∂x³}

(17)

_{u(x - Δx, t) ≈ u(x, t) - Δx ∂u/∂x + 1/2 (Δx)² ∂²u/∂x² - 1/6 (Δx)³ ∂³u/∂x³}

(18)

_{u(x - 2Δx, t) ≈ u(x, t) - 2Δx ∂u/∂x + 2(Δx)² ∂²u/∂x² - 4/3 (Δx)³ ∂³u/∂x³}

(19)

B. Solve KdV Function

% Craig Clark 8/12/18

% Function to solve the KdV equation numerically by an iterative formula.

% This iterative formula is slightly different for the first iteration due

% to the general equation needing the point before as well. It outputs a

% matrix of all of the values of the solution at each grid point at each

% time that it needs to be evaluated.

% Inputs:

% * u_init: The analytical solution is used as the initial values.

% * x: The set of grid points that the solution is to be evaluated at.

% * t: The set of times that the solution is to be evaluated at.

% * delta: The dispersion constant of the system.

% Outputs:

% * u_hist: A matrix of all of the values of the solution at each time.

% Limitations and Constraints:

% - The function will run without any issues, no matter the input but the

% solution will become unstable depending on the input parameters as

% discussed in the script that runs this function.

function [u_hist] = SolveKdV(u_init, x, t, delta)

M = size(t, 2); % Computes the number of times that will be iterated

N = size(x, 2); % Computes the number of grid points

u_hist = zeros(N, M); % Creates a matrix to store all of the values

u_hist(:, 1) = u_init(1:N, 1); % Populates the first column with the initial solution

delta_x = x(2) - x(1); % Finds what the distance step is

delta_t = t(2) - t(1); % Finds what the time step is

u_float = zeros(N + 4, 3); % Creates a float variable to store the necessary

values for the iteration

u_float(:, 1) = u_init(:, 1); % Populates the first column with the initial solution

% Finds the first iteration based on the modified formula

for j = 3:N

u_float(j, 2) = u_float(j, 1) - (delta_t / delta_x) * (1/6) * (u_float(j + 1, 1) +

u_float(j, 1) + u_float(j - 1, 1)) * (u_float(j + 1, 1) - u_float(j - 1, 1)) - delta^2 * (delta_t / (2 * delta_x^3)) * (u_float(j + 2, 1) - 2 * u_float(j + 1, 1) + 2 * u_float(j - 1, 1) - u_float(j - 2, 1));

end

% Uses the periodic boundary conditions to fill in the grid points that cannot

be otherwise calculated

u_float(1, 2) = u_float(N + 1, 2);

u_float(2, 2) = u_float(N + 2, 2);

u_float(N + 3, 2) = u_float(3, 2);

u_float(N + 4, 2) = u_float(4, 2);

u_hist(:, 2) = u_float(1:N, 2); % Attributes these values to the correct part of

the historic values matrix

t1 = 1; % Sets up indexing that can rotate for the float matrix

t2 = 2;

t3 = 3;

for k = 2:M - 1 % Iterates through the main formula for the remaining times to be

calculated

for i = 3:N - 2 % Runs the formula for the current time for each grid point

u_float(i, t3) = u_float(i, t1) - (delta_t / delta_x) * (1/3) * (u_float(i + 1, t2) +

u_float(i, t2) + u_float(i - 1, t2)) * (u_float(i + 1, t2) - u_float(i - 1, t2)) - delta^2 * (delta_t / (delta_x^3)) * (u_float(i + 2, t2) - 2 * u_float(i + 1, t2) + 2 * u_float(i - 1, t2) - u_float(i - 2, t2));

end

% Uses the periodic boundary conditions again

u_float(1, t3) = u_float(N + 1, t3);

u_float(2, t3) = u_float(N + 2, t3);

u_float(N + 3, t3) = u_float(3, t3);

u_float(N + 4, t3) = u_float(4, t3);

u_hist(:, k) = u_float(1:N, t3); % Assigns these values to the correct time

% Cycles through the indexing so that the next time can be calculated

temp = t1;

t1 = t3;

t3 = t2;

t2 = temp;

end

C. Script Animate Solution KdV

% Craig Clark 8/12/18

% Script to solve the KdV equation numerically using the SolveKdV

% function, allowing the user to input the parameters for the time intervals

% at which the solution is computed, how many grid points the x-axis is

% split into and the value of the dispersion constant delta. The script

% will animate the wave, producing a graph at each time step that pauses

% for half a second.

% Constraints and Limitations:

% - This script is limited by the solution becoming unstable as you progress

% to later times.

% - It should also be noted that the higher the resolution the wave is

% evaluated to (the higher the number of grid points), the sooner the

% solution becomes unstable.

% - The value of delta needs to be below around 0.05 otherwise the solution

% becomes unstable almost immediately.

deltat = input('At what time intervals would you like the wave evaluated? ');

N = input('How many grid points would you like the intervals split into? ');

delta = input('What is the value of the dispersion constant, delta? ');

a = delta * sqrt(12); % Sets the scaling for the width of the peak

x = linspace(0, 1, N); % Creates the grid points in the displacement

t = linspace(0, 0.2, 0.2 / deltat + 1); % Finds all of the times at which a solution is computed

M = size(t, 2); % The total number of time steps that will be computed

u_init = zeros(1, N + 4); % Creates an extended vector for all of the initial values

of the wave at each grid point

u_init(1, 1:N) = (sech((x - 0.5) / a)).^2; % Computes the initial values based on the

analytical solution

% Applying the periodic boundary conditions:

u_init(1, N + 1) = u_init(1, 1);

u_init(1, N + 2) = u_init(1, 2);

u_init(1, N + 3) = u_init(1, 3);

u_init(1, N + 4) = u_init(1, 4);

u_init = u_init'; % Turning the vector into a column rather than a row vector

[u_hist] = SolveKdV(u_init, x, t, delta); % Calls up the function to output the value

of u at every time step

figure % Opens up a figure that will animate the wave

grid on % Turns the grid on so that you can see the wave move more easily and make

an estimate of the speed

% The animation of the wave

for n = 1:M

td(1, 1) = deltat * (n - 1); % Calculates the current time for the graph title

titlename = ['Solution to the KdV Equation at t=', num2str(td(1, 1)), 's']; % The

string for the graph title

plot(x, u_hist(:, n)) % Plots the graph at the current time step

grid minor % Turns on the minor grid lines

axis([0 1 -0.5 1.5]) % Fixes the axis so that they don't jump around when the

wave moves/changes amplitude

xlabel('Displacement, x')

ylabel('Diturbance Amplitude, u')

title(titlename)

pause(0.5) % Pauses the loop for half a second so that you can see each step in

time

end

References

1. Oxford University Physics Department CO23 Lab Script: Solitions https://www-teaching.physics.ox.ac.uk/practical course/scripts/srv/local/rscripts/trunk/Computing/CO23/CO23.pdf
2. N.J Zabusky and M.D Kruskal Interaction of ”Solitions” in a Collisionless Plasma and the Recurrence of Initial States https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.15.240
3. Klaus Brauer, The University of Osnabrück, Germany The Korteweg-de Vries Equation: History, exact Solutions, and graphical Representation http://people.seas.harvard.edu/ jones/solitons/pdf/025.pdf