[Arjen Markus] (4 september 2019) Some time ago I came across a nice article about the http://mathworld.wolfram.com/Korteweg-deVriesEquation.html%|%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 https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.15.240%|%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
}
======

<<categories>>Mathematics|Physics