# studydtqconv: Study DTQ Convergence In Rdtq: Density Tracking by Quadrature

## Description

`studydtqconv` facilitates the generation of convergence plots, i.e., plots where one studies the error (in various norms) as a function of the time step h; the error is computed as the difference between the exact PDF and the approximate PDF computed via the DTQ method.

## Usage

 ```1 2``` ```studydtqconv(method, drift, diffusion, exact, hseq, kseq, Mseq, init, fT, thresh = 0) ```

## Arguments

 `method` This must be a string, either "cpp" or "sparse", that indicates which algorithm to use. See the parameter of the same name in the `rdtq` function. `drift` if `method="cpp"`, then this should be a pointer to a drift function implemented in C++ using `Rcpp`. If `method="sparse"`, then this should be the name of an R function. For further details, see the description of the `drift` parameter for the `rdtq` function. `diffusion` if `method="cpp"`, then this should be a pointer to a diffusion function implemented in C++ using `Rcpp`. If `method="sparse"`, then this should be the name of an R function. For further details, see the description of the `diffusion` parameter for the `rdtq` function. `exact` an R function that accepts two arguments, `xvec` and `T`, and returns the exact probability density function p(x,T) for each x in `xvec`. `hseq` a numeric vector of values of h, the time step, to use for the computation of the DTQ solution. Note that `hseq` and `kseq` must have the same lengths. `kseq` a numeric vector of values of k, the grid spacing, to use for the computation of the DTQ solution. Note that `hseq` and `kseq` must have the same lengths. `Mseq` a numeric vector of integer values of M. For each corresponding value of k, the spatial grid will cover the domain [-y_M, y_M] where y_M = Mk. This corresponds to the parameter `bigm` in the `rdtq` function. Note that `kseq` and `Mseq` must have the same lengths. `init` a scalar initial condition. `fT` a positive numeric scalar giving the final time at which to compare the exact and DTQ solutions. `thresh` This optional positive scalar is only used when `method="cpp"`. See the parameter of the same name in the `rdtq` function.

## Value

The function returns the errors between the DTQ and exact solutions indexed by the corresponding value of hseq. The errors are returned in the L^1 norm, L^∞ norm, and the Kolmogorov-Smirnov norm. The errors are returned in the form of a data frame.

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74``` ```# In this example, we will study the convergence of the DTQ method # for the SDE with drift f(x) = x/2 + (1 + x^2)^(1/2) and # diffusion g(x) = (1 + x^2)^(1/2). library(Rdtq) library(Rcpp) # implement the drift and diffusion functions using C++ sourceCpp(code = '#include using namespace Rcpp; double drift(double& x) { return(0.5*x + sqrt(x*x + 1.0)); } double diff(double& x) { return(sqrt(x*x + 1.0)); } typedef double (*funcPtr)(double& x); // [[Rcpp::export]] XPtr driftXPtr() { return(XPtr(new funcPtr(&drift))); } // [[Rcpp::export]] XPtr diffXPtr() { return(XPtr(new funcPtr(&diff))); }') # implement the drift and diffusion functions using R mydrift = function(y) { return(0.5*y + sqrt(y^2 + 1)) } mydiff = function(y) { return(sqrt(y^2 + 1)) } # implement the exact solution at time t, i.e., # the analytical formula for the pdf p(x,t) exactsol = function(xvec,t) { transx = asinh(xvec) - t prefac = (1 + xvec^2)^(-1/2) z = prefac*dnorm(x=transx) return(z) } # define the sequence of parameters that will be used to study convergence hseq = c(0.5,0.2,0.1,0.05) kseq = hseq^(0.55) Mseq = ceiling(5*(-log(hseq))/kseq) # we will use the method="sparse" code for the three largest values in hseq, # and then switch to the method="cpp" code for the three smallest values firstpart = c(1:2) errpart1 = studydtqconv(method="sparse",drift=mydrift,diffusion=mydiff,exact=exactsol, hseq[firstpart],kseq[firstpart],Mseq[firstpart], init=0,fT=1) errpart2 = studydtqconv(method="cpp",drift=driftXPtr(),diffusion=diffXPtr(),exact=exactsol, hseq[-firstpart],kseq[-firstpart],Mseq[-firstpart], init=0,fT=1,thresh=1e-16) # now we will put everything together into one data frame mydat = rbind(errpart1,errpart2) # we plot the convergence diagram, on a log-log scale, using ggplot2 library(ggplot2) library(scales) myplot = ggplot(data=mydat, aes(x=x,y=y,group=norm,color=norm)) myplot = myplot + theme_bw() + theme(plot.background = element_rect(fill='white')) myxticks = sort(10^(round(log(hseq)/log(10)*10)/10)) rawyticks = round(log(mydat\$y)/log(10)*10)/10 rawyticks = round(seq(from=min(rawyticks),to=max(rawyticks),length.out=length(myxticks))*1)/1 myyticks = unique(10^rawyticks) myplot = myplot + scale_x_log10(breaks = hseq) myplot = myplot + theme(axis.text.x = element_text(angle=90,hjust=1)) myplot = myplot + scale_y_log10(breaks = myyticks, labels = trans_format("log10", math_format(10^.x))) myplot = myplot + labs(x="h (temporal step size)", y="error") myplot = myplot + geom_line() + geom_point() # save the plot to a pdf (portable document format) file ggsave(filename="example.pdf", plot=myplot, width=5, height=4) ```