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 }