Residual terms in numerical integration of ODEs - normal error or reality?

[sethoflagos]

I set up an excel spreadsheet to integrate the following system by first order forward difference equations:

 

2/(k-1)*dc/dt + du/dt = 0 in direction dr/dt = u + c

2/(k-1)*dc/dt - du/dt = 0 in direction dr/dt = u - c

 

With r=0 boundary conditions u=Asin(wt), c = c0, I got a solution that was a very good fit to

 

u=Asin(wt)cos(wr/c0) which is the expected result for a standing wave in a cylinder.

 

How good a fit?

 

Well the integration yielded maximum variance values of 0.5 in both r and t trends over one wavelength, and subtracting the above expression reduced this by 99.9954% in the r trends, 99.99902% in the t trends. A pretty good fit.

 

However, looking at the residual and experimenting a little, I found that adding a further term

 

wr(A/c0)^2/pi*sin(2wt)cos(2wr/c0)

 

removed 99.75% of the r, and 98.97% of the t remaining residual max variances.

 

This term may have a very small value over one wavelength, but the proportionality to r means that eventually the second harmonic will dominate the waveform, and the proportionality with A^2 means that it's proportional significance increases with input amplitude.

 

Any suggestions as to how I might go about determining whether this term is no more than an artifact of non-linearities in the forward difference method, or if it really is a true component of the system?

 

In passing, the small residual that remains is almost equal to a simple product of sin(2wt)sin(2wr/c0), leaving something closely proportional to r^2sin(3wt)cos(3wr/c0). And still, little sign of randomicity.

 

Many thanks in advance for your time,

Seth

Edited April 3, 2017 by sethoflagos