Problem simulating a wave on a string by solving the wave equation as a difference equation.

Question

This is a question I posted earlier however the purpose of what I'm doing is going to be expounded much more clearly.

yy[x_, t_] := Which[x == 0, yy[0, t] = 0, x == 1, yy[1, t] = 0, t == 0, yy[x, 0] = Exp[-1000 (x - .3)^2]]

This code is the boundary condition and initial condition for my wave. At x=0 the amplitude equals zero and at x=1 the amplitude is also zero for all times. For t=0 the wave is a Gaussian centered at x=.3

Flatten[Table[yy[x, t + 3./10000.] = yy[x + 1./100., t] + yy[x - 1./100., t], {t,  0, 3./10000., 3./10000.}, {x, 1./100., 99./100., 1./100.}], 1]

This runs my problem and uses the values of yy at t=0 to calculate t=3./10000. later. However when I run this program I get these nulls for the next time step of t=3./10000. This runis the calculation and I cannot get enough good data point for later times.

Help will be very much appreciated. I never had this problem before and I tried everything I know and none of it worked to fix this problem.

Answer 1

It may be overkill, but here's a version with everything rounded that gets rid of the extraneous Null.

Clear[yy];
dx = 1./100; dt = 1./10000;
yy[x0_, t0_] := With[{x = Round[x0, dx], t = Round[t0, dt]},
  Which[
   x == 0, yy[0, t] = 0,
   x == 1, yy[1, t] = 0,
   t == 0, yy[x, 0] = Exp[-1000 (x - .3)^2]]
  ]

Table[With[{
   x = Round[x, dx],
   t = Round[t, dt],
   t2 = Round[t + 3./10000., dt]},
  yy[x, t2] = yy[Round[x + 1./100., dx], t] + yy[Round[x - 1./100., dx], t]
  ],
 {t, 0, 3./10000., 3./10000.},
 {x, 1./100., 99./100., 1./100.}]

See Memoization of Rounded inputs for related strategies.

---

Some explanation

The extraneous values of Null come from Which, not because Which is misused by not having a default case, but because of pattern-matching problems in the memoization of yy. (The OP seems to have noticed something like this in How to stop Mathematica from turning .02 into 0.019999999999999997`, although it is not remarked upon explicitly in the question.)

The underlying cause is round-off error. Note, for instance, that 0.01 + 0.01 equals 0.02`, while 0.03 - 0.01 equals 0.019999999999999997`. This is due to numerical round-off error consistent with IEEE 754 binary64 floating-point numbers. Note that 0.02 - (0.03 - 0.01) equals 2^(-58), a one-bit loss of precision.

The problem with the code is subtler, though. While 0.019999999999999997` == 0.02` returns True, as patterns, 0.019999999999999997` does not match 0.02`. The problem arises specifically in the snippet of the OP's code,

yy[x - 1./100., t]

when x = 0.03. For t = 3./10000, the value has been memoized as yy[0.02`, 0.0003`]; but code snippet evaluates first as yy[0.019999999999999997`, 0.0003`], which in turn evaluates to Null because it does not match yy[0.02`, 0.0003`] and it falls through the Which cases. The rounding in the solution I gave, some of which is not strictly necessary, makes the arguments in the memoized definitions and evaluations of yy match each other when they are meant to be equal.