PDF
The KdV Equation and Solitons
Introduction
The Korteweg-de Vries (KdV) equation, formulated in 1895 by Korteweg and de Vries, models water waves. It contrasts sharply to the Burgers equation because it introduces no dissipation and the waves travel seemingly forever. In 1965, Zabusky and Kruskal named such waves solitons.
The KdV equation with boundary conditions and initial value is formulated as
The equation models the steepening and dispersion of wavefronts but does not support a train of simple harmonic waves. Such trains comprise the wave crests normally associated with the ocean: simply a momentary constructive interference of contributing waves moving at different speeds. However, the equation does support solitons, single “humps” that travel without changing shape or speed for unexpectedly long distances.
Indeed, Perry and Schimke (Ref. 2) concluded from shipboard oceanographic measurements that bands of choppy water in the Andaman Sea, which lies east of the Bay of Bengal and west of Burma (Union of Myanmar) and Thailand, are associated with large-amplitude oceanic internal waves. Satellite images have since clarified that these waves originate on shallow banks on a layer between warm and cool water. Further, Osborne and Burch (Ref. 1) analyzed oceanographic data in an effort to assess the forces of underwater current fluctuations associated with such waves on offshore drilling rigs. They concluded that the visually observed roughness bands are caused by internal solitons that follow the KdV equation (Ref. 3).
A more recent development is the application of the KdV equation to another type of waves — light waves. Today solitons have their primary practical application in optical fibers. Specifically, a fiber’s linear dispersion properties level out a wave while the nonlinear properties give a focusing effect. The result is a very stable, long-lived pulse (Ref. 3). It is amazing that researchers have discovered a formula for such waves:
This equation says that the pulse speed is what determines the pulse amplitude and the pulse width. The following simulation illustrates this effect. An initial pulse, which does not conform to the formula, immediately breaks down into two pulses of different amplitudes and speeds. The two new pulses follow the formula and thus can travel forever. While the formula does not reveal how solitons interact, the simulation shows that they can collide and reappear, seemingly unchanged, just as linear waves do, another counterintuitive observation that is difficult to observe without predictions by computing.
Model Definition
In the model, the term uux describes the focusing of a wave and uxxx refers to its dispersion. The balancing of these two terms permits waves to travel with their shape unchanged.
Because COMSOL Multiphysics does not evaluate third derivatives directly, you rewrite the original equation above as a system of two variables to solve it:
Using the General Form PDE interface, you need to define two dependent variables, u1 and u2, and identify the da, Γ, and F coefficients in the following equation:
•
Only the first equation has a time derivative, and it is with respect to u1, so only da(1,1) is 1; the other three components are zero.
•
The divergence is a space derivative with respect to x. This means that the Γ component from the first equation is u2, which you type as u2. The Γ component from the second equation is u1x, which you express using COMSOL Multiphysics syntax as u1x.
•
The F term components are the right-hand side of the equations: F1 is 6u1u1x (type 6*u1*u1x), and F2 is u2 (type u2).
The initial condition for u1 uses a hyperbolic cosine function to provide an interesting waveform to start with. For u2, you must provide the second space derivative of this function to provide consistent initial conditions.
The boundary conditions are periodic boundary conditions: the solution at one end is always identical to the one at the other end of the domain.
Results
The following plot shows how solitons collide and reappear with their shape intact.
Figure 1: Solution visualizing a soliton collision.
References
1. A.R. Osborne and T.L. Burch, “Internal Solitons in the Andaman Sea”, Science, vol. 208, no. 4443, pp. 451–460, 1980.
2. R.B. Perry and G.R. Schimke, “Large-Amplitude Internal Waves Observed Off the Northwest Coast of Sumatra”, J. Geophys. Res., vol. 70, no. 10, pp. 2319–2324, 1965.
3. G. Strang, Applied Mathematics, Wellesley-Cambridge, 1986.
Application Library path: COMSOL_Multiphysics/Equation_Based/kdv_equation
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  1D.
2
In the Select Physics tree, select Mathematics>PDE Interfaces>General Form PDE (g).
3
Click Add.
4
In the Number of dependent variables text field, type 2.
5
Click  Study.
6
In the Select Study tree, select General Studies>Time Dependent.
7
Click  Done.
Root
1
In the Model Builder window, click the root node.
2
In the root node’s Settings window, locate the Unit System section.
3
From the Unit system list, choose None.
Keeping track of units is not important in this model; by turning off unit support, you avoid the need to specify dimensions for equation coefficients and coordinates to get rid of unit warnings.
Geometry 1
Interval 1 (i1)
1
In the Model Builder window, under Component 1 (comp1) right-click Geometry 1 and choose Interval.
2
In the Settings window for Interval, locate the Interval section.
3
In the table, enter the following settings:
Coordinates
-8
8
General Form PDE (g)
Periodic Condition 1
1
In the Model Builder window, under Component 1 (comp1) right-click General Form PDE (g) and choose Periodic Condition.
2
In the Settings window for Periodic Condition, locate the Boundary Selection section.
3
From the Selection list, choose All boundaries.
General Form PDE 1
1
In the Model Builder window, click General Form PDE 1.
2
In the Settings window for General Form PDE, locate the Conservative Flux section.
3
In the Γ text-field array, type u2 on the first row.
4
In the Γ text-field array, type u1x on the second row.
5
Locate the Source Term section. In the f text-field array, type 6*u1*u1x on the first row.
6
In the f text-field array, type u2 on the second row.
7
Locate the Damping or Mass Coefficient section. In the da text-field array, type 0 in the second column of the second row.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the u1 text field, type -6*sech(x)^2.
4
In the u2 text field, type -24*sech(x)^2*tanh(x)^2+12*sech(x)^2*(1-tanh(x)^2).
Mesh 1
Edge 1
In the Mesh toolbar, click  Edge.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type 0.1.
5
Click  Build All.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.0025,2).
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 3e-6.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
4
From the Method list, choose Generalized alpha.
The Generalized alpha time stepper is well suited for wave problems. For an accurate solution, use tighter tolerance settings.
5
Click to expand the Absolute Tolerance section. From the Tolerance method list, choose Manual.
6
In the Absolute tolerance text field, type 3e-7.
7
Locate the Time Stepping section. In the Amplification for high frequency text field, type 0.98.
8
In the Model Builder window, expand the Study 1>Solver Configurations>Solution 1 (sol1)>Time-Dependent Solver 1 node, then click Fully Coupled 1.
9
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
10
In the Tolerance factor text field, type 0.1.
11
In the Study toolbar, click  Compute.
Results
1D Plot Group 1
1
In the Settings window for 1D Plot Group, locate the Data section.
2
From the Time selection list, choose From list.
3
In the Times (s) list, select 0.025.
Line Graph 1
1
In the Model Builder window, expand the 1D Plot Group 1 node, then click Line Graph 1.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type -u1.
4
In the 1D Plot Group 1 toolbar, click  Plot.
The solution to the KdV equation at 0.025 s.
To visualize the solution, extrude results along the time axis.
Parametric Extrusion 1D 1
In the Results toolbar, click  More Datasets and choose Parametric Extrusion 1D.
2D Plot Group 2
In the Results toolbar, click  2D Plot Group.
Surface 1
1
Right-click 2D Plot Group 2 and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type -u1.
Height Expression 1
1
Right-click Surface 1 and choose Height Expression.
2
Click the  Zoom Extents button in the Graphics toolbar.
Compare with the plot shown in Figure 1.