Message Boards Message Boards

1
|
2315 Views
|
9 Replies
|
5 Total Likes
View groups...
Share
Share this post:

How can I best mimic Fortran Double Precision in MMA ?

I am trying to reproduce some results that were originally generated by a Fortran package that executed on a CDC 7600 mainframe (back in the late 1970s). The results are obtained by summing a large number (up to 15000) rapidly oscillating terms in what is called a Partial Wave Sum.

I cannot get my MMA code to match the old results and I think that it may be down to numerical rounding/convergence issues. I believe that Fortran double precision meant that a DP variable was stored as two physical words of memory (each of which had 60 bits on a 7600).

I'm running MMA on a MacBook Pro with 64-bit Apple silicon and I have tried tweaking Working Precision, Accuracy Goal and the like but with no visible benefits.

Is it possible to coerce MMA into using two 64-bit words per variable, or are there more elegant ways of handling these kinds of issue?

Any help gratefully received.

David

POSTED BY: David Mackay
9 Replies

Hello Eric,

Thanks again for your comments and suggestions, I really appreciate all the time you have spent looking at my examples. I need to take time to experiment with the ideas shared but I will post a response to let you know if I overcome my current issues.

Hans - thank for the tip on editing the overall thread title.

Bye, David

POSTED BY: David Mackay

Many thanks for the help and tips so far.

I have attached a Notebook which illustrates the basic problem. The key numbers that I am trying to calculate are generated in the module etaLIST and are stored in the list etas. At the end of the Notebook you can see how the precision of these data varies with the parameter l.

In my real problem, I have to sum up a series which may contain up to 15000 terms. Each term has three parts. One part is the exponent of an eta and so the precision of the final sum is critically dependent upon the accuracy and precision with which I can compute the etas themselves.

I may have gone too far with my use of backticks !!

David

Attachments:
POSTED BY: David Mackay
Posted 1 year ago

Just started looking, but I've already seen something I want to comment on. Integers and rationals are generally handled as perfect/infinite precision quantities by Mathematica. So, for example, optEXP = 20`50; is actually less precise than just optEXP = 20;. If you want high precision, one strategy is to stick with infinite precision as long as possible, but here you've already "corrupted" any calculations using optEXP. [Note, obviously you may need to trade speed for precision.]

So, this:

sK = 2.46227`50;
sA = 53.401191`50;
optC = 200.0`50;
optEXP = 20`50;
wK = optC/sA^2;

could be this instead:

sK = Rationalize[2.46227`50];
sA = Rationalize[53.401191`50];
optC = 200;
optEXP = 20;
wK = optC/sA^2;

And this:

voverE[r_] := (4`30/sK)*(1`50/r^12 - 1`50/r^6) - (I*wK)*(1`50/r)^optEXP

could be this:

voverE[r_] := (4/sK)*(1/r^12 - 1/r^6) - (I*wK)*(1/r)^optEXP
POSTED BY: Eric Rimbey
Posted 1 year ago

Another comment (but this isn't really my area of expertise, so someone else may have more insight)... it looks like you've gone to a lot of trouble to define certain quanties with 50 digits of precision, but in your FindRoot calculations, you are setting WorkingPrecision to 25. Furthermore, the documentation states that AccuracyGoal and PrecisionGoal by default are set to WorkingPrecision/2. So, I think that we can only expect a dozen or so digits of precision in the results. Is that what you want?

Also, after waiting for several minutes for the final two calculations to complete, I just gave up. So I don't know what you're trying to show. Please provide a more focused and manageable example of where you're losing precision unexpectedly (sorry, but I'm just not gonna invest my whole day in this).

POSTED BY: Eric Rimbey
Posted 1 year ago

We could better help you if you provided some sample code that represents what you want to do. In the meantime, I'll guess that you are running into a common challenge for new Mathematica users trying to control precision. Whenever you write a number like 1.0, i.e. anytime you use a decimal point and no other precision info, you are specifying a number with machine precision (about 16 decimal digits). Once you've committed to a particular precision, you cannot expect precision to increase in subsequent calculations. Setting the working precision for a particular calculation cannot increase the precision of the quantities you feed to that calculation. So, the general strategies are to use infinite precision quantities for your inputs (things like integers or Pi) or to specify the precision of your inputs sufficiently high. One way to specify precision is with a backtick, e.g. 1.0`100 (approximately 100 digits of precision in this case). But again, we'd need a concrete example of what you're trying to do to provide more specific guidance.

POSTED BY: Eric Rimbey

Hello Eric, many thanks for your responses. I haven't shared a Notebook yet as the results I am attempting to match are, sadly, paper based and so you would not be able to see the problem I have in trying to match them. I can upload a jpeg of the original data though. Before I do that though, I am experimenting with the use of ` to specify desired precisions and starting with 30 which would seem to be double (not noble, sorry about missing the spell corrector blunder) precision on my Mac.

If I try this, I can set a high precision value for an estimate for FindRoot for example. The result returned, however, still seems to be machine precision. I guess that I am missing a step here and I would certainly appreciate further help.

I can build a stripped down .nb to illustrate my attempts if that would help but you may be able to suggest some more things for me to try first?

Many thanks, David

POSTED BY: David Mackay
Posted 1 year ago

Have you read the documentation for FindRoot? There are examples of how to use AccuracyGoal, PrecisionGoal, and WorkingPrecision. I am able to get results to very high precision by tweaking those options. Please don't provide a complicated notebook. Just provide a single, simple, representative example of something that doesn't give you what you expect.

POSTED BY: Eric Rimbey
Posted 1 year ago

Off-topic, but:

(not noble, sorry about missing the spell corrector blunder)

You can correct a typo in the title of a discussion you have started. Hit the Edit button below your published post. Then you can edit the title, in the input field on top.

enter image description here

POSTED BY: Hans Milton

Welcome to Wolfram Community! Please provide your efforts in the form of the Wolfram Language code. This will make it easier for other members to help you. Check several methods available to include your code in the rules http://wolfr.am/READ-1ST

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract