Hello, I am testing a well accepted 6th order compact scheme. The first and second order norm is shown at www.glue.umd.edu/~tgma

Form that curve, we can see that the scheme produce the 6th order accuracy at first 3 points, then the error is almost constant no matter what gridsize you use. It maybe come from the truncation error from the machine, which I used Dell PC 700MHz. The compiler is Digital Fortran, and it is Fortran90 code. I wonder if anybody touched similiar topic, and can provide me the reasons for that phenomena? Thanks in advance! Tingguang

Are you running and printing with double precision ?

Yes, most variable is Real(8). For output, I used f20.15 as the format. The L2 norm stablizes after a certain gridnumber, that is what I am wondering. Is it truncation error? But I have went beyond 1e-8 before, that is, I got L2 norm 1e-11 before. So I cannot explain this result here. Thanks if any one can provide any suggestions! Tingguang

The numerical constants in the scheme should also be double precision.

E.g.:

x = 1.0d0/3.0d0 * y

instead of

x = 1.0/3.0 * y

Regards, Markus

Markus is on the right track - the whole calculation should be done in double precision (real*8) unless you are _really really_ sure the single precision part will _never ever_ limit your accuracy. The problem with the L2 norm staying larger than expected is most likely due to round-off errors, not truncation errors, and the level of these is problem-dependent. If you get L2 down to 1.d-11 in one problem, that is no guarantee it will go that low in the next problem.

Thanks, Markus and Dean, I rechecked the code. Basically, I follow the rules on double procision. But I found that the constant Pi (to use periodic boundary condition, I have to use Pi), I used only 8 digits before, When I change to 23 digits, it produced the following result(second diagram). www.glue.umd.edu/~tgma When I calculated the slope of Log10(L2 norm) vs. Log10(1/delta x), the slope is found to be 5.83 (only for the first 5 points), while my scheme is 6th order accuracy. I wonder if it is the final result, even I have only 5 points agree with prediction? When is the minimum error I can have (for this PC)? Can anyone provide any insight for this roundoff error? Thanks in advance! Tingguang

This does look like a typical precision problem. Maybe your machine does not really give double precision. I know that for example my Sun workstation's double precision is not at all as good as the NEC super-computer that I have access to. It would be a good idea to try it in another machine.

How are you calculating the error? Do you have an analytical solution? Since you have already got 5.83rd order accuracy, it is surely some external factor that is affecting it.

chidu...

The empirical 5.83 order may be as good as it gets. The formal order of a method is determined in most cases by a Taylor series expansion of the difference equations. That formal order is an upper bound on the actual behavior under grid refinement for a variety of reasons. First, the total error is the sum of truncation and round-off errors. The former get smaller as the grid is refined, the latter get worse. At some level of resolution, the round-off becomes larger and the solution stops improving with further grid refinement. Second, even if you could ignore the round-off, that formal order is the actual order only in the limit of grid size going to zero. If one considers what multiplies the delta x^6 in the context of the Taylor series with remainder, it is typically a sixth derivative or a product of high-order derivatives evaluated somewhere on the grid, and that "somewhere" changes with delta x. Finally, if the solution has discontinuities or steep gradients (shocks, flames, and the like) or solutions with singularities in the lower order derivatives, the empirical order of the method can be lowered way below 6. As a simple example, try using Simpson's rule to integrate x**(n/2) from 0 to 1 with n=1, 3, 5. Evaluate the empirical order of the integrations, and you will see n=1 produces a reduction of at least one order. The problem is due to the singularity at 0 of the first derivative. This is a subtlety that seems to be ignored in most numerical analysis courses and texts.

Although you are right to question the validity of your result and search for an error or other explanation, don't be surprised if 5.83 is as good as it gets for this one test problem.

For periodic problems with a finite range of length scales, I would not be surprised to find an exact 6th order behaviour given enough precision.

chidu...

Thanks, chidu and Dean, I have another 8th order scheme test result posted @ www.wam.umd.edu/~tgma . The calculated precision is 7.57 (for only first 4 points). It seems that high order scheme is good where total error level is larger than 1e-8. Go beyond that limit, the improvement is very limited. This high order compact scheme is said to have an accuracy comparable with Spectral-like resolution, anyone had any idea toward that? My test fuction is only sin and cos, very simple and periodic.

Tingguang

Chidu,

I agree that with the idealizations you list, you might get the full 6th order. My point was that in the real world, the range of length scales is determined by the physics (and may be effectively infinite), and that can degrade the order. Also, real world computers don't let one use infinite precision in the arithmetic operations, and round-off will always be there. The question is whether or not it is significant in any particular simulation. This is something we can't determine for the case that began this thread; the person who did that calculation is going to have to do numerical experiments with the same code, computer, and test problem to find out whether the results are affected by round-off.

There are modules in f90 available that define new class types to give you the precision you want. I think the basic premise is that for quad precision you just store the variable as two separate ones and use operator overloading to add, subtract, etc..

you just need the space to store the variables

hmm...sorry I think this is off topic

You are absolutely right! Now I can get my error norm to 1e-14. Thanks! Tingguang

I'm afraid that I can't help you with your specific problem, but there is a nice article entitled What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg which explains round-off errors etc. quite well - might help you to track down the source of your problems.