Numerical Glitch - could someone please advise

Hello Guys, I've been working on a rather complicated code for quite some time. Recently I noticed some rather strange behaviour in my results, however, and have been running tests to see if I can identify where the problem is.

I can't post the code here, which sadly means I can't really talk about the exact details of what the code is doing. I think I can explain the problem, however, without going into specifics - i'd really appreciate it if someone could help me work out where I'm going wrong.

Basically the code calculates a large array of numbers, the case I'm dealing with should result in a perfectly symmetrical array; the algorithm requires two steps, the first step is an initial approximation (which should be symmetrical) and the second step is a refinement of the first. The problem is that I am noticing asymmetry within the results - this seems to grow with time and is very worrying. I've specifically designed the code to produce exactly symmetric results in each step and I am getting very small errors which eventually get rather large.

The problems always seem to occur in the same region of my array, near the top and bottom edges. To test for symmetry I subtract the top value in a column from the bottom element in that column. If the values cancel then my results are good - if not then something is going wrong. In fortran my code looks like

Step 1 - phi1 -> phi2

Test Symmetry Via
   DO I=1,imax
         IF (abs(phi2(I,jmax)-phi2(I,1))>0) THEN
           write(*,*) 'problem found, timestep ', step, ' in column I=', I,'
         ENDIF 
   ENDDO

Step 2 phi2 -> phi3

where phi1 is the initial array (the old value), phi2 the intermediate and phi3 is the final array. The code uses a timestepping procedure and the equations are constructed so that the array should be exactly symmetrical (essentially it does the same set of calculations twice for each column and assigns the values symmetrically down each column of the array).

The problem is that if I test the symmetry between step 1 and 2 then the act of testing seems to change my output. if I run the code, letting it do step 1, then step 2 without the testing procedure then I get a different output than if I include the procedure. Commenting out the IF statement seems to change the output data.

The test is an exact quote from my code - it's just an IF statement that checks if two values are equal or not and doesn't change any variables. I don't know why this should have an effect on the output.

If someone could please help me I'd greatly appreciate it. This problem is a bit weird and is driving me a tad crazy 😮)

Many thanks, B[/code]

A couple of ideas that might apply.

I'd check how phi1 and phi2 are defined, so there is no memory overwrite that can occur if they are passed as subroutine arguments.

Also check the range of definition for phi2, in case you are not defining all the values.

If a finite difference itteration, is the edge rule being correctly applied.

Otherwise write out the phi2 updates to a text file and import them into a spreadsheet to check. ( Excel 2007 can support some fairly large data sizes now !)

best of luck

John C

If the calculation is sensitive to round-off error then you can get different results when IO statements etc are inserted. This is because the change can affect the optimisation process and hence change the round-off error.

FTN95 has compile time switches that allow you to change all REAL4 variables to REAL8 or even REAL*10 and this could provide some insight into the effect of round-off error.

However, if your calculation is sensitive to round-off error then you will probably not be able to put any confidence in the results and you may have to investigate the numerical stability of the algorithm that you are using.

Surely this statement is very sensitive to precision:-

IF (abs(phi2(I,jmax)-phi2(I,1))>0) THEN 

as the slightest difference in either of the two variable values can swap the logic between false and true.

Instead of comparing against zero, wouldn't a very small value be better, say:-

abs(phi2(I,jmax)+phi2(I,1))/1.0d8

or whatever you judge to be a small enough number ?

Maybe but I was thinking of something else.

If you have a number of assignment statements and you interpose an output statement (that should in theory have no effect on the assigned values), then the assigned values can sometimes change because of a change in the underlying compiler optimisations. The compiler optimised computations can have different round-off errors. Hence, if what you are getting depends only on the round-off error then you will get different results.

This may sound strange but this kind of effect is not unusual in very long calculations and in short calculations that are numerically unstable.

Thanks guys, I've made sure that phi1, and phi2 are filled matrices, the errors I'm finding are very small 1e-16 etc. and i'm pretty sure they're just rounding errors.

My problem is that these errors are asymmetrically distributed down each column. I know that this seems a tad strict but I have tried to force this by performing each calculation twice and allocating the results to each column. Any asymmetry in the results seems to grow quite rapidly.

In my case I have an array of 250 columns and 60 rows. In my tests I have forced the first 100 columns and final 50 columns to be symmetric by assignments like

phi1(1:30,1:100)=phi1(60:31:-1,1:100) 

and so on. I'm only solving for the middle columns - this seems to be the region in which I'm encountering difficulties.

The algorithm requires each column to be solved sequentially, each column requires 2 passes with a tridiagonal solver. One pass solves phi(1:30,I), the second provides phi(60:31:-1,I), i.e. from the outer edges towards the centre. The top and bottom halves of the array are uncoupled, hence solving separately. The equations are constructed to ensure symmetric distribution of rounding errors - I expect a completely symmetric solution.

I have checked the edge conditions and ensured that my equations are diagonally dominant.

I realise that it does seem a bit strict to want perfect symmetry but after extensive testing I've realised that the problem is extremely sensitive and my results are not good without it. Solving half the domain, then reflecting the results, is not an option here.

Going crazy.... B