Large numerical errors in my Mathematica program

Question

Here is a piece of my program. Have a look.

For[m = 1, m <= mode1, m++,
  For[n = 0, n <= mode2, n++,
    A[m, n][t_] = a[m, n]*Cos[\[Omega]*t];
    B[m, n][t_] = b[m, n]*Cos[\[Omega]*t];
  ]
]

temp = 0;
For[m = 1, m <= mode1, m++,
  For[n = 0, n <= mode2, n++,
    temp++;
    equation[temp] = 
      ExpandAll[Integrate[eqC[m, n]*Cos[\[Omega]*t], {t, 0, (2*Pi)/\[Omega]}]];
    equation[temp] = ExpandAll[Simplify[equation[temp]/10^9]];
    Print["\n\nEquation ", temp, "-\n", equation[temp]];
    temp++;
    equation[temp] = 
      ExpandAll[Integrate[eqS[m, n]*Cos[\[Omega]*t], {t, 0, (2*Pi)/\[Omega]}]];
    equation[temp] = ExpandAll[Simplify[equation[temp]/10^9]];
    Print["\n\nEquation ", temp, "-\n", equation[temp]];
  ]
]

After running of this code I am supposed to get few equations and then create a matrix out of it by a series of differentiations. I know that the matrix must come out to be symmetric. The problem is that when I enter simple data i.e. e=1,h=1, etc. I get accurate results and the matrix is symmetric, but as soon as I give the real data which have values like 71.02e9,0.000247 the calculations come out to be wrong and I get an unsymmetric matrix. I have thoroughly checked the code and cannot find a single mistake on my part. I have also checked the results of the program for a simple case with manual calculations.

Your post doesn't really contain a clear question, nor have you provided enough information for us to run your code (which is not written in the normal Mathematica style - it's too imperative). However, your problem is probably just loss of [precision](http://reference.wolfram.com/mathematica/guide/NumericalEvaluationAndPrecision.html). The two "real data" values you gave will be interpreted as `MachinePrecision` numbers. Either use exact numbers or bump the precision up to say 30 by using ``71.02`30*^9`` and ``0.000247`30``. You should also try to use the built-in Fourier transform functions — Simon, Dec 16 '11 at 05:43
So you've got all these As and Bs at the start but not in the next block. And you refer to e=1 and h=1 even though e and h aren't anywhere in the code. At a guess it looks like you are trying to solve for Fourier series coefficients or something??? — Mike Honeychurch, Dec 16 '11 at 05:51
Yes I agree that the code that I have written might not be the most elegant way to program in mathematica. I am actually new to the language. I am not trying to develop a Fourier transform, I am solving a partial differential equation of shell vibration using galerkin's method. Also the 'e' and the 'h' are actually in the eqC[m,n] and eqS[m,n] which have been declared earlier. — , Dec 16 '11 at 12:01
Let me rephrase what Mike asked, what variables in the above code are you setting to your "real values?" It is not clear which ones they are. — rcollyer, Dec 16 '11 at 12:20
The `ExpandAll[Simplify` step just after the previous `ExpandAll` seems superfluous to me. You can add the division by 10^9 directly to the first step. — Sjoerd C. de Vries, Dec 16 '11 at 12:35

Timo · Answer 1 · 2011-12-16T11:33:04.533

7

You can try and increase the precision of your calculations by globally setting, e.g., $MinPrecision=50 and specifying your data values to a high precision using either foo = SetPrecision[0.000247,50] or using the shorthand 0.000247`50.

edited Dec 16 '11 at 11:33

answered Dec 16 '11 at 05:42

Timo

4,246
6
29
42

score 2 · Answer 2 · answered Dec 16 '11 at 08:18

As a further alternative you could use Rationalize[0.000247] and then derive a numerical quantity later with N[expr, prec]. In M- if you give inexact input you get inexact output and for exact input you get exact output.

Large numerical errors in my Mathematica program

2 Answers2