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

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.

• I provided a brief explanation of the problem and a link to a solution in comments to your other question. Did you look at them? – Michael E2 Apr 7 '15 at 0:43
• I did see them. I don't see how rounding is going to help me. My numbers are less then far less then one so when I round them they just turn to zero. – Daniel Berkowitz Apr 7 '15 at 0:48
• Did you try Round[t, 0.0001] or something like that? – Michael E2 Apr 7 '15 at 0:50
• You mean yy[Round[x,t]]? – Daniel Berkowitz Apr 7 '15 at 0:55
• No, I meant to round x and t to multiples of 0.01 and 0.0001, which seems appropriate. That would take care of the round-off error. But trying it out, doesn't seem to work. There may be something about your scheme I don't understand yet. – Michael E2 Apr 7 '15 at 0:58

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.

• It works thank you so much I can't express how happy I am. I'll study your code intensively so that this problem does not happen again. – Daniel Berkowitz Apr 7 '15 at 1:08
• @DanielBerkowitz You're welcome. – Michael E2 Apr 7 '15 at 1:09
• Personally, I feel this Q&A should be merged with the OP's other one, 79130. – Michael E2 Apr 7 '15 at 2:44
• @DanielBerkowitz FWIW, the same computation can be achieved with nsteps = 2; NestList[Join[{0.}, ListConvolve[{1, 0, 1}, #], {0.}] &, DeveloperToPackedArray@Table[N@yy[x, 0.], {x, 0., 1., 1./100.}], nsteps] - It's about ten times faster or more. – Michael E2 Apr 7 '15 at 12:48