10
$\begingroup$

I am quite a newbie and it is taking me so much time to understand what I need to do to achieve something really simple and I still haven't found it.

I run a dynamic model of a metabolic pathway in Mathematica and I calculate the values of several variables at steady-state for multiple conditions. In my model I have parameters that are constant. These have low precision just because it would cause quite some effort to put a bunch of zeros after each one (and it also makes the code look unappealing). (There must be a way I think to tell Mathematica to consider these numbers as high precision numbers, right?)

So for all the different conditions I calculate the variable values at steady-state. And I use ParallelTable quite often and I also use FindRoot and NDSolve.

An example of a relevant piece of script is:

tsol0 = ParallelTable[
   NDSolve[Join[
     Odes /. RateEqs /. CoAMATX /. 
       ReplacePart[Parm, 
        Position[Parm, C6AcylCarMAT][[1, 1]] -> 
         C6AcylCarMAT -> k[[j]]] /. Parm2, InitialConditions], 
    Vars, {t, 0, 1000000}], {j, 1, Length[k]}];

ssp0 = Table[
   FindRoot[
    Table[Odes[[i, 2]] == 0, {i, 1, Length[Odes]}] /. RateEqs /. 
       CoAMATX /. 
      ReplacePart[Parm, 
       Position[Parm, C6AcylCarMAT][[1, 1]] -> 
        C6AcylCarMAT -> k[[j]]] /. Parm2, 
    ParallelTable[{Vars[[i]][t], (Vars[[i]][900000] /. tsol0[[j]])[[
       1]]}, {i, 1, Length[Vars]}]], {j, 1, Length[k]}];

Parm and Parm2 contain lists of constant parameters. Odes, RateEqs and CoAMATX are functions and equations that together define how all the variables should change in time. With tsol I calculate the values for all the variables for each time point. I do this for k conditions, in each of these conditions the constant parameter C6AcylCarMAT has different value, because I want to calculate the variable values in time for an increasing C6AcylCarMAT value. With ssp0 I want to find the steady-state value for each of the conditions by indicating that probably at around = 900000 mins the variable values should not change anymore.

As a result, ssp0 is a table containing the steady-state values for all the (metabolite concentration-related) variables for each condition. (Each row is one specific condition and each column is one specific metabolite). The other steady-state variables I can calculate from the steady-state metabolite concentrations and Odes, RateEqs and CoAMATX.

Anyway long story. This all works, however I get precision-related warnings. Furthermore, my output values do not have the same amount of decimals. This is the warning/error I get below my ssp0 code:

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

Plus I believe that I have huge error propagations, especially because some of the steady-state variable values are very large or very small and thus further calculations with these values cause quite some errors. And when I do my control analysis calculations, only in the conditions where some of the variable values are very large or very small, the checks don't come out well; so some of the values the identity matrix that need to be zero are not, they're like 0.6 etc instead of zero which is a huge error.

I just want to set at the beginning of the notebook:

  1. That I want to set all the constant parameters (in Parm and Parm2 e.g.) at a very high precision (for instance 30)
  2. That I want all the calculations in the notebook to be done with high precision
  3. That I want all the results to be displayed with a high precision (thus e.g. 30 decimals).

Can you guys please help me? I have been trying lots of things and breaking my head for quite some time and I gave up now. I also tried 

$PreRead = (# /. 
     s_String /; 
       StringMatchQ[s, NumberString] && 
        Precision@ToExpression@s == MachinePrecision :> 
      s <> "`50." &);

from the thread "How to set the working precision globally? $MinPrecision does not work" but that also didn't work (at least it did not display the results with all these decimals, so now i don't know for sure whether it worked or not).

Improve this question
$\endgroup$
7
  • $\begingroup$ If you would please try for format your question as stipulated here it would be great :) $\endgroup$
    – Öskå
    Oct 6, 2014 at 16:06
  • $\begingroup$ It is really a good idea to use symbol names starting with a lowercase letter not uppercase. i.e. RateEqs bad; rateEqs good. All of Mathematica's own symbols start with uppercase, if you try to use one of them as your own symbol you'll run into some serious problems. $\endgroup$
    – Ymareth
    Oct 6, 2014 at 16:11
  • $\begingroup$ Dear Öskå, Thank you for your suggestion. Introduction: In my notebook I define and analyse a dynamic model of a biochemical pathway: Each reaction, uses a unique set of parameters (listed in parm and parm2). The reaction rate equations are listed in rateEqs. I use odes to define how the metabolite concentrations change in time. I define the variable metabolites (vars) and their concentration at time zero (initialConditions). I then use tsol and ssp to calculate concentrations in time and at steady-state for k different substrate concentrations (C6-CarMAT). $\endgroup$
    – Anne-Claire
    Oct 8, 2014 at 7:47
  • $\begingroup$ Dear Ymareth, Thank you for your suggestion, I will change my symbols. $\endgroup$
    – Anne-Claire
    Oct 8, 2014 at 7:49
  • $\begingroup$ Analysis Method in the notebook: For the model analysis I also calculate tsol and ssp for the k substrate concentrations, for different parameter values. For instance, I vary all the parm2 with lets say 1%. With these results I calculate a variable called control coefficient (CC). In addition I calculate a variable called elasticity (Elast) by calculating the partial derivatives of the rate equations. Each of these variables I calculate for all the k substrate concentrations. Then I use CC and Elast to calculate matrices called Ci and E. $\endgroup$
    – Anne-Claire
    Oct 8, 2014 at 7:59

2 Answers 2

Reset to default
8
$\begingroup$

You can increase the minimum precision used in a notebook by setting $MinPrecision to a value greater than 0 (the default). For instance, if you want to ensure that all numbers are defined to have precision twice that of machine precision, you could specify:

$MinPrecision = 2 $MachinePrecision;

To make sure that your parallel kernels also use this same precision, run the above command in a ParallelEvaluate call, before you run any other parallel commands:

ParallelEvaluate[$MinPrecision = 2 $MachinePrecision];

I find this is the easiest way to guarantee higher precision is used in most calculations.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Dear Guillochon, Thank you for this suggestion. Can I put this at the top of my script and it will work for all my ParallelTable calculations? $\endgroup$
    – Anne-Claire
    Oct 8, 2014 at 8:24
  • $\begingroup$ @Anne-Claire It should, so long as you don't later kill the parallel kernels. The ParallelEvaluate needs to be run every time the parallel kernels are created. $\endgroup$
    – Guillochon
    Oct 9, 2014 at 0:39
7
$\begingroup$

Real numbers generally have $MachinePrecision

x = .25 // Precision

MachinePrecision

You can designate arbitrary precision for a number using a number mark and designated precision

x = .25`10; y = .25`30;

Precision /@ {x, y}

{10., 30.}

You can get a rational approximation with Rationalize

.25 // Rationalize

1/4

However, Rationalize doesn't always work unless you specify a tolerance. Use tolerance of 0 to force a rational output

x = 2.897654321*^-16 // Rationalize

2.89765*10^-16

x = 2.897654321*^-16 // Rationalize[#, 0] &

1/3451067274494265

In general, Rationalize[#, 0]& all of your real or complex constants except for exact symbolic constants (e.g., Pi, E). Rationals are exact numbers and have infinite precision. Any subsequent numerical operations can be set to whatever required WorkingPrecision and the constants will not cause precision warnings.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ RootApproximant can be considered as an alternative to Rationalize, but unless it is really expected inputs to be either Rationals or Algebraics, these are a hacks. Arbitrary precision arithmetic correctly tracks precision of results on basis of precision of input (it typically reduces over computations), unlike computations with "precise" inputs. $\endgroup$
    – kirma
    Oct 7, 2014 at 5:48
  • $\begingroup$ Dear Bob, Thank you for these tips. I have used `100 and also WorkingPrecision and SetPrecision to increase the precision of the parameters and to increase the working precision for the tsol and ssp calculations. However, my problem still exists. None of my results changed, other than that they have more decimals. $\endgroup$
    – Anne-Claire
    Oct 8, 2014 at 8:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.