need help on testing scheme and error
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 
Re: need help on testing scheme and error
Are you running and printing with double precision ?

Re: need help on testing scheme and error
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 1e8 before, that is, I got L2 norm 1e11 before. So I cannot explain this result here. Thanks if any one can provide any suggestions! Tingguang

Re: need help on testing scheme and error
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 
Re: need help on testing scheme and error
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 roundoff errors, not truncation errors, and the level of these is problemdependent. If you get L2 down to 1.d11 in one problem, that is no guarantee it will go that low in the next problem.

Re: need comments!
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

Re: need comments!
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 supercomputer 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... 
Re: need comments!
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 roundoff errors. The former get smaller as the grid is refined, the latter get worse. At some level of resolution, the roundoff becomes larger and the solution stops improving with further grid refinement. Second, even if you could ignore the roundoff, 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 highorder 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. 
Re: need comments!
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... 
More comments!
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 1e8. Go beyond that limit, the improvement is very limited. This high order compact scheme is said to have an accuracy comparable with Spectrallike resolution, anyone had any idea toward that? My test fuction is only sin and cos, very simple and periodic.
Tingguang 
Re: need comments!
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 roundoff 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 roundoff. 
Re: comment!
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 :( 
Re: need help on testing scheme and error
You are absolutely right! Now I can get my error norm to 1e14. Thanks! Tingguang

Re: need comments!
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 FloatingPoint Arithmetic by David Goldberg which explains roundoff errors etc. quite well  might help you to track down the source of your problems.

All times are GMT 4. The time now is 10:48. 