Update README.md

README.md

@@ -12,16 +12,66 @@ elliptic, convection-diffusion type of partial-differential equation
 
 ```math 
 	\forall x_{\text{\tiny H}}\in \mathbb{R}^n: \begin{cases}
-	x_{\text{\tiny H}}^\intercal Qx_{\text{\tiny H}} - J + x_{\text{\tiny H}}^\intercal A^\intercal \nabla V + \frac{1}{2}\text{Tr}(R\nabla^2V) = 0, \\  
-	V(x_{\text{\tiny H}})\leq \rho + V(0), 
+	x_{\text{\tiny H}}^\intercal Qx_{\text{\tiny H}} - J + x_{\text{\tiny H}}^\intercal A^\intercal \nabla V + \frac{1}{2}\text{Tr}(R\nabla^2V) = \frac{\partial V}{\partial t}, \\  
+	V(x_{\text{\tiny H}},t)\leq 0, 
 	\end{cases}\quad\quad
 	 \forall x_{\text{\tiny H}} \in \partial \Omega: 
 	\begin{cases}
-	V(x_{\text{\tiny H}}) = \rho + V(0),\\
+	V(x_{\text{\tiny H}},t) = 0,\\
 	\nabla V = 0.
 	\end{cases}
 ```
 
-The solution to this PDE is the value function $`V`$, and the free boundary
-$`\partial \Omega`$ is the threshold on the state $`x_{\text{\tiny H}}`$ which 
-defines the optimal sampling scheme. For more details, we refer to the [article](cdc2017_paper.pdf).
+The solution to this PDE is the value function $`V`$. In stationarity, the free boundary
+$`\partial \Omega`$ is the threshold on the vector state $`x_{\text{\tiny H}}`$ which 
+defines the optimal sampling scheme. The current implementation handles second-order systems.
\ No newline at end of file
+For more details, we refer to our [article](cdc2017_paper.pdf).
+
+# Running the Code
+For getting started, we recommend running and modifying the script `script_optimal_bound_2D.m`. The stationary solution of $`V`$ is computed by running the function `optimal_bound_2D`, found in `functions/optimal_bound_2D.m`. A typical run is:
+
+```matlab
+%Define system
+A = [-1 20; 20  -1];
+q11 = 1;
+q22 = 1;
+J=1;
+
+%Gridding
+dx = 3e-2;
+x1min=-2;
+x1max = 2;
+x2min= -2;
+x2max = 2;
+x1vec = x1min:dx:x1max;
+x2vec = x2min:dx:x2max;
+h = 5e-4;
+
+%Enable/disable plotting and printing
+plot_progress = 0;
+print_progress = 1;
+
+%Compute V
+V = optimal_bound_2d(A, q11, q22, J, x1vec, x2vec, h, plot_progress, print_progress);
+```
+Without loss of generality, `optimal_bound_2D` assumes that the system is of a form where $`R = I`$ and $`Q = \text{diag}(q_{11}, q_{22})`$ (the system state can always be transformed to achieve this).
+The input arguments of `optimal_bound_2D` are:
+
+* `A` : The $`2\times2`$ system matrix of the linear reset system.
+* `q11`,`q22` : The diagonal elements of the (diagonal) cost matrix Q.
+* `J` : The total cost of state deviation and sampling.
+* `x1vec`,`x2vec` : The grid vectors of the two state dimensions. Each vector is assumed to have equidistantly spaced elements.
+* `h` : The time step of the simulation.
+* `plot_progress` : Specifies if $`V`$ and $`\partial \Omega`$ should be plotted during the simulation (`plot_progress = 1`) or not (`plot_progress = 0`).
+* `print_progress` : Specifies if iteration info should be printed to the command window during the simulation (`print_progress = 1`) or not (`print_progress = 0`).
+
+# About the Numerical Method
+The numerical solution method is based on simulating the non-stationary version of the PDE using finite-difference approximations of the differential operators and gridding in time and space. 
+In this implementation, we use the backward-time central-space (BTCS) finite-difference approximation (see e.g. this [article](http://depa.fquim.unam.mx/amyd/archivero/DiferenciasFinitas4_25730.pdf) for more details).
+
+We simulate the PDE by iterating the BTCS scheme, while at each iteration enforcing the inequality $`V \leq 0`$ by the assignment `V = min(V,0)`. This has the effect of moving the
+free boundary $`\partial \Omega`$ over time. The simulation is then progressed until a stopping criterion based on stationarity of $`V`$ is met, or if the maximum number of iterations
+is reached. The stopping criterion is an upper threshold on the maximum absolute difference of $`V`$ between to successive iterations.
+
+The current implementation is for second-order systems only, i.e. $`x_{\text{\tiny H}} \in \mathbb{R}^2`$. However, this method of applying BTCS can be directly generalized to higher order systems
+as well, although at a higher computational cost (the time complexity scales exponentially with the state dimension). 
\ No newline at end of file