## Solitons and the Korteweg-de Vries equation

Arjen Markus (4 september 2019) Some time ago I came across a nice article about the Korteweg-de Vries equation , an equation arises in various physical contexts which has the remarkable property of being both non-linear and allowing for waves that do not lose their identity - so-called solitons. The article actually describes experiments with a variant, but that does not quite matter. What does matter is that it had a receipe for solving the equation numerically, so that I did not have invent one myself ;).

The program below simply implements the receipe and prints three time series. While experimenting with it, I found that the method requires very small time steps, because at some point a parasitical solution seems to take over and the calculation explodes. A small time step reduces the growth of that solution. I have not implemented a graphical display of the resulting curve yet, but that should be a fairly straightforward exercise. Although, given the small time step requirment and the subsequent fairly long calculation time, it also provides an opportunity to demonstrate the use of multithreading - one thread doing the calculation, the other presenting the result in a graph.

In the meantime, havea look at the article.

```# kdv.tcl --
#     Small program inspired by a paper by Zabusky and Kruskal (https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.15.240)
#     to solve the Korteweg-de Vries equation:
#                   2
#     u  + u u   + d  u    = 0
#      t      x        xxx
#
#     The parameter d can chosen freely, but the program follows the paper: d = 0.022
#     Also:
#     - The interval is 0 <= x <= 2
#     - The initial condition is a cosine function with period 2
#     - The boundary conditions are cyclic
#

# nextSolution --
#     Calculate the solution for the next time step
#
# Arguments:
#     unow           The current solution
#     uprev          The previous solution
#
# Result:
#     The solution for the next time step
#
# Note:
#     All parameters are set as global variables
#
proc nextSolution {unow uprev} {
global d
global dx
global dt

set n [llength \$unow]

set advec [expr {\$dt / \$dx / 3.0}]
set diff  [expr {\$d**2 * \$dt / \$dx**3}]

set unext {}
for {set i 0} {\$i < \$n} {incr i} {
set upr [lindex \$uprev \$i]

set im2 [expr { (\$i-2) % \$n }]
set im1 [expr { (\$i-1) % \$n }]
set ic  [expr { (\$i)   % \$n }]
set ip1 [expr { (\$i+1) % \$n }]
set ip2 [expr { (\$i+2) % \$n }]

set um2 [lindex \$unow \$im2]
set um1 [lindex \$unow \$im1]
set uc  [lindex \$unow \$ic ]
set up1 [lindex \$unow \$ip1]
set up2 [lindex \$unow \$ip2]

set du [expr {\$upr - \$advec * (\$um1 + \$uc + \$up1) * (\$up1 - \$um1)
- \$diff  * (\$up2 - 2.0*\$up1 + 2.0*\$um1 - \$um2) }]

lappend unext \$du
}

return \$unext
}

# main ---
#     Set up the parameters and run the program for a while
#
set n 100
set d 0.022

set pi [expr {acos(-1.0)}]
set dx [expr {1.0 / \$n}]
set dt 1.0e-5

set unow {}
for {set i 0} {\$i < 2*\$n} {incr i} {
lappend unow [expr {cos(\$pi*\$i*\$dx)}]
}

set uprev \$unow

for {set t 0} {\$t < 2000000} {incr t} {
set unext [nextSolution \$unow \$uprev]

if { \$t % 10000 == 0 } {
puts "\$t\t[lindex \$unext [expr {\$n/4}]]\t[lindex \$unext [expr {\$n/2}]]\t[lindex \$unext [expr {3*\$n/4}]]"
}

set uprev \$unow
set unow  \$unext
}```

 Category Mathematics Category Physics