両端を固定した弦のシミュレーションになっています。黒い線が現在の状態とすると、薄い線は一つ前のステップです。

【使い方】

1. 弦は２０区間に差分しています。自由な代表点の数は１９です。初期設定は本来自由なのですが、今回、そこには凝っていません。で、コンボボックスの３つのメニューから選べます。
2. 「three waves」というのは、３つの波頭という意味ではなく、フーリエ成分の３つを適当に合成して作ったものです。
3. 「cembalo」というのは、弦を一点で引っ張った初期状態から始めたもので、チェンバロやギターの状況に近いと考えます。
4. 「piano」は逆に、初期の変位はないが、一定範囲に初速度を与えた条件になっています。フェルトのハンマーを使った実際の状況に近いかどうかは保証の限りでないです。 明らかな違いは、こちらは本当に減衰が少ないということです（勿論、これから減衰を加えることはできます）。
5. 「single」を押すと、時間のワンステップだけ動き、「count」が１つ増えます。
6. ４０ステップで波が一巡し、最初の状態に戻ります。ただし、最初の状態を記憶して戻しているわけではなく、実際にシミュレーションしています。 「multi100」を押すと、4千ステップ、すなわち、波の100巡分だけ計算して表示します。変化がないように見えると思います。
「multi1e5」を押すと、4百万ステップ、すなわち、波の10万巡分だけ計算して表示します。これも、変化がないように見えます。
7. 「to-left」をチェックすると、今回の代表点と前回のすぐ右側の代表点とを繋ぐ線分が表示されます。 「single」を押して変化を見ると、変化分が右側、つまり、川上（最初「川下」としていましたが、間違い。「下から上へ」、です）に移動しているのがわかると思います。位置の上下ではなく、傾きに着目してください。
8. 「to-right」をチェックすると、今回の代表点と前回のすぐ左側の代表点とを繋ぐ線分が表示されます。 「single」を押して変化を見ると、変化分（傾き）が左側（やはり、川上）に移動しているのがわかると思います。
9. 「anime」をチェックすると、動画のように動きます。

以上の理論です。

一次元の波動方程式は次の形を持つ。
1 s² ∂²u ∂ t² = ∂²u ∂ x² 　　　　（1）
シミュレーションにおいては、時間および空間を離散化して扱うが、それぞれの間隔をt_Pおよびx_Pとすると、
T ≡ t t_P 　　　　（2）
X ≡ x x_P 　　　　（3）
ここに T および X は、シミュレーションにおいて、通常は整数とされる数値である。さらに、
S ≡ t_P s x_P 　　　　（4）
と置くと、
1 S² ∂²u ∂ T² = ∂²u ∂ X² 　　　　（5）
時間に対する微係数は、離散化すると、
∂²u ∂ T² ≒ 1 2 (u_T=i+1,X － 2 u_T=i,X + u_T=i-1,X) 　　　（6）
この式は、時刻依存性が放物線であることが保証される場合に限れば正確に成り立つ。一方、時刻依存性が正弦曲線で、サンプル時点間の位相がθ_Tである場合は、次のようにすべきである。
∂²u ∂ T² f(θ_T) = 1 2 (u_T=i+1,X － 2 u_T=i,X + u_T=i-1,X) 　　　（7）
ここに、
f(θ) = 1 - cos(θ) θ²　　　（8）
は補正係数であり、せいぜい f(θ_T) < π の範囲で有効である。
また、位置座標に対する微係数も同様に離散化することができ、
∂²u ∂ X² ≒ 1 2 (u_T,X=j+1 － 2 u_T,X=j + u_T,X=j-1) 　　　（9）
これも、正確なのは、位置依存性が放物線になる場合のみである。仮に曲線が正弦曲線で、サンプル位置間の位相がθ_Xである場合は、次のようにすべきである。
∂²u ∂ T² f(θ_X) = 1 2 (u_T=i+1,X － 2 u_T=i,X + u_T=i-1,X) 　　　（10）
ここに、f(θ) は、先に（8）として定義された補正係数である。これらの置換えで（5）式は、
u_T=i+1,X=j = 2 u_T=i,X=j － u_T=i-1,X=j + S² f(θ_X) f(θ_T) (u_T=i,X=j+1 － 2 u_T=i,X=j + u_T=i,X=j-1)　　　（11）
ここに、θ_X や θ_T は、対象となる波の波長や振動数が明らかでないと決まらないものなので、 そもそも色々な波長が混ざった現実の波動をシミュレートするする手段としては不十分であるように見える。
一つの解決は、空間分割、時間分割ともに、解析する波長、振動数よりも十分に短く取ることである。このとき、 f(θ_X)、f(θ_T) ともに１に近づくので、その比もほぼ１になる。 注意点として、S² は１よりは大きくない値にする必要があって、さもなくば、時間を経るにつれ、誤差による振動が発生する。
もう一つの解決は、S² を、１、または、１に極めて近い１未満の値にすることである。 波動は、波長によらず、一定速度で伝わるので、空間のアトムと時間のアトムとの比をこの速度（位相速度）にすることで、 θ_X と θ_T とを常に等しくできる。 そのことで、荒いメッシュであっても、ぎりぎりの振動数まで精度の高いシミュレーションを提供する。
空間のメッシュの取り方には、例えば弦の長さの整数分の１にするなどの制約があろうから、離散化間隔は、時間について最終調整するのが現実的だろう。
仮に、S² f(θ_X) f(θ_T) を１とした場合、上の式はこうなる。
u_T=i+1,X=j = u_T=i,X=j+1 + u_T=i,X=j-1 － u_T=i-1,X=j　　　（12）
（5）式で X = 0, N を固定端とし、 S = 1 とした場合の解として、次のフーリエ級数があり、
u(T, X) = A sin(π m T/N) sin(π m X/N) 　　　（13）
一般解はその m=1～∞ の和であると言えるが、m < N の範囲の（13）式は、(12)式を満たす。
実際、（13）を（12）に代入するなら、左辺は、
A sin((i+1) m π/N) sin(j m π/N)　　　（14）
であって、右辺は、
A sin(i m π/N) sin((j+1) m π/N) + A sin(i m π/N) sin((j-1) m π/N) － A sin((i-1) m π/N) sin(j m π/N)　　　（15）
　　= A [sin(i m π/N) {sin(j m π/N) cos( π/N) + sin(m π/N) cos(j m π/N)}
　　　+ sin(i π/N) {sin(j m π/N) cos(m π/N) － sin(m π/N) cos(j m π/N) }
　　　－ {sin(i m π/N) cos(m π/N) － sin(m π/N) cos(i m π/N)} sin(j m π T/N) ]
= A { sin(i m π/N) sin(j m π/N) cos(m π/N) + sin(m π/N) cos(i m π/N) sin(j m π/N)}
= A { sin(i m π/N) cos(m π/N) + sin(m π/N) cos(i m π/N)} sin(j m π/N)
= A sin((i+1) m π/N) sin(j m π/N)　　　（16）
となり、確かに左辺に等しいことが示されるのである。
(11) 式から、sin(i m π/N) sin(j m π/N) という中心の項を早々に消去してしまっていたが、ここで示されているのは、個々の代表点を基準とした、X, T 両方向の「曲率」が常に等しいという事実である。
N = 30 で、m = 1, 2, 3 の各次数のフーリエ成分について、３次元グラフを示すが、このことを感じて頂けるだろうか？

　　　

Sや分割数を変えてどの位の精度が得られるか、実際にシミュレーションして計測してみた。
下の表で S=1 としたのが、今回の方法であり、地を赤く塗っている。分割数を上げてこれに匹敵する精度を得るには、利用する最も短い波の１００分の１まで分割しなくてはらならない（青地）。 計算の安定性を考えると、時間の分割もそれ以上にしなくてはならず、このことから、今回の考慮によれば、そうでない場合の少なくとも１万倍の計算速度が得られることがわかる。

表：波長による伝搬速度の変化