Mathematica implementation of the Relaxation algorithm
======================================================
(by Manuel Bichsel, ETH Zurich)
Relaxation algorithm
--------------------
A detailed description of the Relaxation algorithm to numerically solve boundary value problems can be found in 'Press W., Flannery B., Teukolsky S., Vetterling W. (1989): Numerical Recipes in Pascal (Cambridge University Press), p. 645' or in the paper by T. Trimborn, K.-J. Koch and T. M. Steger mentioned below. In short, the algorithm works by starting with an almost arbitrary initial trajectory in phase space and iteratively adjusting this path until it is as close to the solution path as wished. In each iteration the algorithm discretizes time along the momentary trajectory in phase space and computes the real difference between consecutive points on the momentary trajectory. It then also computes the predicted difference between consecutive points on the momentary trajectory as it should be based on the differential equation system. The algorithm then computes the error, i.e. the difference between the real and the predicted difference, which depends on the momentary trajectory. The relaxation algorithm then constructs with Newton's root finding algorithm leading to a linear equation system a new trajectory, whose error in real and predicted differences is (hopefully) smaller than before. These iterations continue until the error gets small enough, as specified by you.
Mathematica implementation
--------------------------
These Mathematica implementations of the Relaxation algorithm for three different economic models (Ramsey-Cass-Koopmans, Jones, Lucas) are based on corresponding implementations done in Matlab by Timo Trimborn, Karl-Josef Koch and Thomas M. Steger (2004), see http://www1.uni-hamburg.de/IWK/trimborn/relaxate.htm, where an accompanying paper is also available.
Caution:
1. The equilibrium based on the ODEs must be finite.
2. The Mathematica implementation can not deal with negative or zero values in the solution, i.e. all values must be positive at all times!
Each of the three Mathematica implementations is available in a version for Mathematica 5 and Mathermatica 6 and consists of four parts: a short description, an input, a computation and an output part. The computation and output parts of all three implementations are exactly the same for each Mathematica version and usually need not be changed. Users should only need to adapt the input part to their model. But users should still check the computation and output part for error and warning messages after each program run.
In the input part you can define parameters and calculate auxilliary variables. This part consists of three sections: model settings, algorithm settings and plot settings.
model settings:
---------------
The following variables and functions _must_ be provided:
- variableNames (list variable): contains the names of all variables. The list begins with "t", and is followed by variable names in the same order as their time derivatives appear in fODE[] (see below).
- fODE[var_1_,var_2_,...,var_n_] (function list with n arguments): contains the right-hand sides of the ODEs, i.e. the time derivatives of the variables in variableNames (without "t") in the same order. The number of arguments n equals nInitial+nFinal+nStatic (see below). The variable list of arguments var_1_, var_2_, ..., var_n_ must also be filled in.
- nInitial, nFinal, nStatic (numerical variables): contain the number of initial and final boundary values, respectively, and the number of static equations (see fSTAT below).
- initialValues and finalValues (list variables): contain the initial and final boundary values in the same order as in variableNames (without "t"). The first nInitial (see above) values of initialValues will be initial boundary values, and the nFinal values of finalValues following directly after the first nInitial ones will be final boundary values (but see also finalEqIndex below). finalValues will also be needed to construct an initial trajectory in phase space which will then be iteratively refined to get the definitve solution.
The following variables and functions _can_ be provided:
- fSTAT[var_1_,var_2_,...,var_n_] (function list with n arguments): contains the static equations, which must be equal to zero at all times. The variable list of arguments var_1_, var_2_, ..., var_n_ must also be filled in.
- finalEqIndex (list variable): final boundary conditions need not be given directly by values but may instead be given by the requirement that some of the ODEs shall be equal to zero at the last time point. In this case the variable finalEqIndex lists all these ODEs by giving their number in the function list fODE, where 1 is the first function, i.e. the time derivative of the first variable after "t" in variableNames. The number of ODEs designated by finalEqIndex must equal nFinal (see above). Caution: finalValues has still to be provided, as it will be needed to construct an initial trajectory in phase space which will then be iteratively refined to get the definitve solution.
algorithm settings:
-------------------
The following variables and functions _must_ be provided:
- timeFunc[tau_] (function): contains the time function regulating the flow of time. Time functions must be invertable and differentiable and map the interval [0,1] to [0,infinity] or [0,tEnd], where 0