Solving continuous models with dependent uncertainty: A computational approach
J.-C. Cortes, J.-V. Romero, M.-D. Roselló, F.-J. Santonja and R.-J. Villanueva
November 27th, 2013


In this website, we are going to present the source Mathematica code [1] and the examples in the paper Solving continuous models with dependent uncertainty: A computational approach, co-authored by J.C. Cortés, J.V. Romero, M.D. Roselló, F.J. Santonja and R.J. Villanueva published in the Journal Abstract and Applied Analysis DOI: 10.1155/2013/983839.

From this link, you can download the main Mathematica code with its structure. Unpack it and it will appear the following structure, under the folder "MAIN".



The main file is the Mathematica notebook "DependentPC.nb". In the folder "Code" the are some text files containing the Mathematica code. Short explanations of the code, step by step, are included. In the folder "Model" we have to include the files describing the model and its equations. In the folder "Output" the procedure will write the output files.

Now, we present, step by step, how to prepare and run the 1st example of the paper. We have to create two files "ModelData.m" and "Data.m" where the characteristics of the model will be described.

The file “ModelData.m” for the 1st example contains the equation of the model and the variables. In general, this file can be also created using epiModel [2], but this 1st example is not long.

(* File ModelData.m *)
(* Equations *)
eqns = { X'[t] + a X[t]^2 + b X[t] + c == 0,
         X[t0] == X0 };

(* Variables *)
Fvars = { X[t] };
vars = { X };

Note that 'a', 'b', 'c' and 'X0' are parameters that we have to determine in the file "Data.m", as constants or random variables following a certain distribution.

(* File Data.m *)
(****************************************************************************************************************
    Variables to generate symbolically the auxiliary system of differential equations.
 ****************************************************************************************************************)
(* Order of the model, i.e., highest derivative of the model *)
ordEDO = 1;
(* Random polynomials truncation order *)
nTrunc = 10;
(* Random variables, grouped by their mutual dependence *)
Grupo = { { b, X0 }
        };

(****************************************************************************************************************
    Distributions of the random variables.
    There are as much distributions as subgroups (sublists) in the variable 'Grupo'.
 ****************************************************************************************************************)
fD = { PDF[ MultinormalDistribution[ {0, 1}, { {3, -1},
                                               {-1, 2} }],
            { Subscript[\[Xi], 1], Subscript[\[Xi], 2] } ]

     };   

(****************************************************************************************************************
    Value of constant model parameters.
    Initial and final time.
 ****************************************************************************************************************)
OtrosDatosModelo = {
    (* Constant parameters *)
    a -> 0,
    c -> 0,
    (* Initial time *)
    t0 -> 0,
    (* Final time *)
    tFin -> 10
           };     

(*****************************************************************************************************************
    Auxiliary variables. DO NOT TOUCH.
 ****************************************************************************************************************)
(* Dimension of the chaos ... *)
dimCaos = Length[ Flatten[ Grupo ] ];
(* ... and random variables *)
varAl = Table[ Subscript[ \[Xi], i ], { i, dimCaos } ];

The above files are copied to "Model" folder. Then, open with Mathematica the file "DependentPC.nb" and execute the cell in the notebook sections:

1.- Generation of the auxiliary system of differential equations

In this section we generate the auxiliary system of differential equations symbolically and we write it in the file "SEDO" inside the folder "Output". Also, we store the list of symbolic inner products to be carried out in the file "listaProductos" inside the folder "Output". Moreover, we calculate the symbolic mean and standard deviation of the solution of the auxiliary system and store them in the file "Media-DT" inside the folder "Output".
Note that, in this section, we built the solutions symbolically, only manipulating expressions. Any computation has been done yet.

2.- Manipulation of the distributions and generation of the inner products

In this section we generate the inner product using the joint probability distributions. The result is written in the file "Inner_Product" inside the folder "Output". As before, any computation has been done yet.
At this point, we should say that if we could simplify the expression of this inner product, further computations will be faster. Therefore, we recommend the user to take advantage of his/her knowledge of Mathematica in order to find simpler forms of the inner product using the ideas of the Section 3.1 of the paper.

3.- Computation of the inner products.

Here, we start with real computations. Despite these computations are parallelized, this is a crucial point and usually consumes most of the procedure computation time. In this section, we compute the inner products from their symbolic expressions stored in the file "Output/listaProductos". It usually takes very long time, the more time the more dependent random variable together (in the same subgroup) because the inner product has more nested integrals.

4.- Resolution of the auxiliary system of differential equations.

The last step of the procedure is to substitute the results of the inner products into the symbolic auxiliary system of differential equations and solve it. We should take into account that, even though our objective has been to work symbolically as much as possible, here we use the Mathematica command NDSolve[], that solves numerically the differential equations. Furthermore, we also have to realise that the auxiliary system uses to have a large number of equations what is an additional difficulty when computing the solution.

5.- Drawing the solutions.

At this point, the procedure has finished. Our last step is to draw the mean and the standard deviation. Some commands appear in the section to plot the results.

More examples

In the following, links to the "ModelData.m" and "Data.m" files of the remainder paper examples are presented. Remember that these files have to be copied inside the "Model" folder. Then, repeat the execution process using the Mathematica notebook "DependentPC.nb".

REFERENCES

  1. Mathematica, http://www.wolfram.com/mathematica/
  2. epiModel, http://epimodel.imm.upv.es