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
}