# R/GQD.passage.R In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions

```GQD.passage<-function(Xs,B,theta,t,delt)
{
#utils::globalVariables(c('fpt','ieq','solverFPT'),'DiffusionRgqd')
fpt      =function(res,res2,N,delt){}
ieq      =function(res1,res2,N,delt){}
solverFPT=function(Xs,Bt,theta,N,delt,tt){}
rm(list =c('fpt','ieq','solverFPT'))
#solverFPT=function(Xs,Bt,theta,N,delt,tt){}
#library(Rcpp)
txt1=
'

arma::vec f(arma::vec a,arma::vec theta,double t)
{

arma::vec atemp(4);

atemp(0)=(+(theta[1])*1+(theta[2])*a(0)+(theta[3])*(a(1)+a(0)*a(0)));
atemp(1)=(+(theta[2])*2*a(1)+(theta[3])*(2*a(2)+4*a(0)*a(1))+(theta[4])*1+(theta[5])*a(0)+(theta[6])*(a(1)+a(0)*a(0)));
atemp(2)=(+(theta[2])*3*a(2)+(theta[3])*(3*a(3)+6*a(0)*a(2)+6*a(1)*a(1))+(theta[5])*3*a(1)+(theta[6])*(3*a(2)+6*a(0)*a(1)));
atemp(3)=(+(theta[2])*4*a(3)+(theta[3])*(8*a(0)*a(3)+24*a(1)*a(2))+(theta[5])*6*a(2)+(theta[6])*(6*a(3)+12*a(0)*a(2)+12*a(1)*a(1)));
return atemp;
}

// [[Rcpp::export]]
{

double val=exp(-0.5*log(2*3.141592653589793*xx(1))-0.5*(Bt-xx(0))*(Bt-xx(0))/xx(1));
return(val);
}
// [[Rcpp::export]]
{

double val=-exp(-0.5*log(2*3.141592653589793*xx(1))-0.5*(Bt-xx(0))*(Bt-xx(0))/xx(1))*(Bt-xx(0))/xx(1);
return(val);
}

// [[Rcpp::export]]
double pcurr(arma::vec xx,double Bt,arma::vec theta)
{
return(val);
}

// [[Rcpp::export]]
arma::vec solverFPT(double Xs,double Bt,arma::vec theta,int N,double delt,double tt)
{
arma::vec fx0(4);
arma::vec fx1(4);
arma::vec fx2(4);
arma::vec fx3(4);
arma::vec fx4(4);
arma::vec fx5(4);
arma::vec fx6(4);
arma::vec fx7(4);
arma::vec fx8(4);
arma::vec fx9(4);
arma::vec fx10(4);
arma::vec fx11(4);
arma::vec fx12(4);
arma::vec fx13(4);
arma::vec fx14(4);
arma::vec fx15(4);
arma::vec fx16(4);
arma::vec x0(4);
arma::vec x1(4);
arma::vec x2(4);
arma::vec x3(4);
arma::vec x4(4);
arma::vec x5(4);
arma::vec x6(4);
arma::vec x7(4);
arma::vec x8(4);
arma::vec x9(4);
arma::vec x10(4);
arma::vec x11(4);
arma::vec x12(4);
arma::vec x13(4);
arma::vec x14(4);
arma::vec x15(4);
arma::vec x16(4);
x0.fill(0);
x0(0)=Xs;
arma::vec fxx0(4);
arma::vec fxx1(4);
arma::vec fxx2(4);
arma::vec fxx3(4);
arma::vec fxx4(4);
arma::vec fxx5(4);
arma::vec fxx6(4);
arma::vec fxx7(4);
arma::vec fxx8(4);
arma::vec fxx9(4);
arma::vec fxx10(4);
arma::vec fxx11(4);
arma::vec fxx12(4);
arma::vec fxx13(4);
arma::vec fxx14(4);
arma::vec fxx15(4);
arma::vec fxx16(4);
arma::vec xx0(4);
arma::vec xx1(4);
arma::vec xx2(4);
arma::vec xx3(4);
arma::vec xx4(4);
arma::vec xx5(4);
arma::vec xx6(4);
arma::vec xx7(4);
arma::vec xx8(4);
arma::vec xx9(4);
arma::vec xx10(4);
arma::vec xx11(4);
arma::vec xx12(4);
arma::vec xx13(4);
arma::vec xx14(4);
arma::vec xx15(4);
arma::vec xx16(4);
xx0.fill(0);
xx0(0)=Bt;
arma::vec res1(N);
arma::vec res2(N);

double valA;
double valB;
double smtemp;
double d=tt;
for (int i = 1; i < N; i++)
{

fx0=f(x0,theta,d);
x1=x0+delt*(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
x10=x0+delt*(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
x11=x0+delt*(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
x12=x0+delt*(0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
x13=x0+delt*(0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15);
fx16=f(x16,theta,d+delt);

x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt;

valA=pcurr(x0,Bt,theta);

fxx0=f(xx0,theta,d);
xx1=xx0+delt*(0.1*fxx0);
fxx1=f(xx1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
xx2=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1);
fxx2=f(xx2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
xx3=xx0+delt*(0.202259190301118*fxx0+0.606777570903354*fxx2);
fxx3=f(xx3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
xx4=xx0+delt*(0.184024714708644*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3);
fxx4=f(xx4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
xx5=xx0+delt*(0.0879007340206681*fxx0+0.410459702520261*fxx3+0.482713753678866*fxx4);
fxx5=f(xx5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
xx6=xx0+delt*(0.085970050490246*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-0.0731856375070851*fxx5);
fxx6=f(xx6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
xx7=xx0+delt*(0.120930449125334*fxx0+0.260124675758296*fxx4+0.0325402621549091*fxx5-0.0595780211817361*fxx6);
fxx7=f(xx7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
xx8=xx0+delt*(0.110854379580391*fxx0-0.0605761488255006*fxx5+0.321763705601778*fxx6+0.510485725608063*fxx7);
fxx8=f(xx8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
xx9=xx0+delt*(0.112054414752879*fxx0-0.144942775902866*fxx5-0.333269719096257*fxx6+0.49926922955688*fxx7+0.509504608929686*fxx8);
fxx9=f(xx9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
xx10=xx0+delt*(0.113976783964186*fxx0-0.0768813364203357*fxx5+0.239527360324391*fxx6+0.397774662368095*fxx7+0.0107558956873607*fxx8-0.327769124164019*fxx9);
fxx10=f(xx10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
xx11=xx0+delt*(0.0798314528280196*fxx0-0.0520329686800603*fxx5-0.0576954146168549*fxx6+0.194781915712104*fxx7+0.145384923188325*fxx8-0.0782942710351671*fxx9-0.114503299361099*fxx10);
fxx11=f(xx11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
xx12=xx0+delt*(0.985115610164857*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-1.37896486574844*fxx5-0.861164195027636*fxx6+5.78428813637537*fxx7+3.28807761985104*fxx8-2.38633905093136*fxx9-3.25479342483644*fxx10-2.16343541686423*fxx11);
fxx12=f(xx12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
xx13=xx0+delt*(0.895080295771633*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3-0.851236239662008*fxx5+0.398320112318533*fxx6+3.63937263181036*fxx7+1.5482287703983*fxx8-2.12221714704054*fxx9-1.58350398545326*fxx10-1.71561608285936*fxx11-0.0244036405750127*fxx12);
fxx13=f(xx13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
xx14=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1+0*fxx2+0*fxx3-0.777333643644968*fxx4+0*fxx5-0.0910895662155176*fxx6+0.0910895662155176*fxx12+0.777333643644968*fxx13);
fxx14=f(xx14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
xx15=xx0+delt*(0.1*fxx0-0.157178665799771*fxx2+0.157178665799771*fxx14);
fxx15=f(xx15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
xx16=xx0+delt*(0.181781300700095*fxx0+0.675*fxx1+0.34275815984719*fxx2+0*fxx3+0.259111214548323*fxx4-0.358278966717952*fxx5-1.04594895940883*fxx6+0.930327845415627*fxx7+1.77950959431708*fxx8+0.1*fxx9-0.282547569539044*fxx10-0.159327350119973*fxx11-0.145515894647002*fxx12-0.259111214548323*fxx13-0.34275815984719*fxx14-0.675*fxx15);
fxx16=f(xx16,theta,d+delt);

xx0=xx0+(0.0333333333333333333333333333333333333333333333333333333333333*fxx0
+0.0250000000000000000000000000000000000000000000000000000000000*fxx1
+0.0333333333333333333333333333333333333333333333333333333333333*fxx2
+0.000000000000000000000000000000000000000000000000000000000000*fxx3
+0.0500000000000000000000000000000000000000000000000000000000000*fxx4
+0.000000000000000000000000000000000000000000000000000000000000*fxx5
+0.0400000000000000000000000000000000000000000000000000000000000*fxx6
+0.000000000000000000000000000000000000000000000000000000000000*fxx7
+0.189237478148923490158306404106012326238162346948625830327194*fxx8
+0.277429188517743176508360262560654340428504319718040836339472*fxx9
+0.277429188517743176508360262560654340428504319718040836339472*fxx10
+0.189237478148923490158306404106012326238162346948625830327194*fxx11
-0.0400000000000000000000000000000000000000000000000000000000000*fxx12
-0.0500000000000000000000000000000000000000000000000000000000000*fxx13
-0.0333333333333333333333333333333333333333333333333333333333333*fxx14
-0.0250000000000000000000000000000000000000000000000000000000000*fxx15
+0.0333333333333333333333333333333333333333333333333333333333333*fxx16)*delt;

d=d+delt;
valB=pcurr(xx0,Bt,theta);

res1(i)=valB;
if(i>1)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+res2(j-1)*res1(i-j);
}
res2(i) = (2.0*valA-2.0*smtemp*delt);
}
}
return(res2);
}

'

txt2=
'

arma::vec f(arma::vec a,arma::vec theta,double t)
{

arma::vec atemp(4);

atemp(0)=(+(theta[1])*1+(theta[2])*a(0)+(theta[3])*(a(1)+a(0)*a(0)));
atemp(1)=(+(theta[2])*2*a(1)+(theta[3])*(2*a(2)+4*a(0)*a(1))+(theta[4])*1+(theta[5])*a(0)+(theta[6])*(a(1)+a(0)*a(0)));
atemp(2)=(+(theta[2])*3*a(2)+(theta[3])*(3*a(3)+6*a(0)*a(2)+6*a(1)*a(1))+(theta[5])*3*a(1)+(theta[6])*(3*a(2)+6*a(0)*a(1)));
atemp(3)=(+(theta[2])*4*a(3)+(theta[3])*(8*a(0)*a(3)+24*a(1)*a(2))+(theta[5])*6*a(2)+(theta[6])*(6*a(3)+12*a(0)*a(2)+12*a(1)*a(1)));
return atemp;
}

// [[Rcpp::export]]
{
double p=(1.0/3.0) *(3*(xx(3)/6.0)*xx(1) - pow(xx(2)/2.0,2))/pow(xx(3)/6.0,2);
double q=(1.0/27.0)*(27*pow(xx(3)/6.0,2)*(xx(0)-Bt) - 9*(xx(3)/6.0)*(xx(2)/2.0)*xx(1) + 2*pow(xx(2)/2.0,3))/pow(xx(3)/6.0,3);
double chk=pow(q,2)/4.0 + pow(p,3)/27.0;
double th=-(xx(2)/2.0)/(3*(xx(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));

double K =xx(0)*th+(xx(1)*th*th)/2.0+(xx(2)*th*th*th)/6.0 +(xx(3)*th*th*th*th)/24.0;
double K1=xx(0)   +(xx(1)*th)       +(xx(2)*th*th)/2.0    +(xx(3)*th*th*th)/6.0;
double K2=xx(1)   +(xx(2)*th)       +(xx(3)*th*th)/2.0;
double val=exp(-0.5*log(2*3.141592653589793*K2)+(K-th*K1));
return(val);
}
// [[Rcpp::export]]
{
double p=(1.0/3.0) *(3*(xx(3)/6.0)*xx(1) - pow(xx(2)/2.0,2))/pow(xx(3)/6.0,2);
double q=(1.0/27.0)*(27*pow(xx(3)/6.0,2)*(xx(0)-Bt) - 9*(xx(3)/6.0)*(xx(2)/2.0)*xx(1) + 2*pow(xx(2)/2.0,3))/pow(xx(3)/6.0,3);
double chk=pow(q,2)/4.0 + pow(p,3)/27.0;
double th=-(xx(2)/2.0)/(3*(xx(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
double thdash=1.0/(xx(1)+th*xx(2)+0.5*th*th*xx(3));

double K =xx(0)*th+(xx(1)*th*th)/2.0+(xx(2)*th*th*th)/6.0 +(xx(3)*th*th*th*th)/24.0;
double K1=xx(0)   +(xx(1)*th)       +(xx(2)*th*th)/2.0    +(xx(3)*th*th*th)/6.0;
double K2=xx(1)   +(xx(2)*th)       +(xx(3)*th*th)/2.0;
double gg      = 1.0/sqrt(2*3.141592653589793*(K2));
double ggdash = -3.141592653589793*pow(2*3.141592653589793*(K2),-3.0/2.0)*(xx(2)*thdash+xx(3)*thdash*th);
double hh      = exp(K-th*K1);
double hhdash = exp(K-th*K1)*(-th*thdash*xx(1)-th*th*thdash*xx(2)-0.5*(th*th*th)*thdash*xx(3));

double val=(ggdash*hh+gg*hhdash);
return(val);
}

// [[Rcpp::export]]
double pcurr(arma::vec xx,double Bt,arma::vec theta)
{
return(val);
}

// [[Rcpp::export]]
arma::vec solverFPT(double Xs,double Bt,arma::vec theta,int N,double delt,double tt)
{
arma::vec fx0(4);
arma::vec fx1(4);
arma::vec fx2(4);
arma::vec fx3(4);
arma::vec fx4(4);
arma::vec fx5(4);
arma::vec fx6(4);
arma::vec fx7(4);
arma::vec fx8(4);
arma::vec fx9(4);
arma::vec fx10(4);
arma::vec fx11(4);
arma::vec fx12(4);
arma::vec fx13(4);
arma::vec fx14(4);
arma::vec fx15(4);
arma::vec fx16(4);
arma::vec x0(4);
arma::vec x1(4);
arma::vec x2(4);
arma::vec x3(4);
arma::vec x4(4);
arma::vec x5(4);
arma::vec x6(4);
arma::vec x7(4);
arma::vec x8(4);
arma::vec x9(4);
arma::vec x10(4);
arma::vec x11(4);
arma::vec x12(4);
arma::vec x13(4);
arma::vec x14(4);
arma::vec x15(4);
arma::vec x16(4);
x0.fill(0);
x0(0)=Xs;
arma::vec fxx0(4);
arma::vec fxx1(4);
arma::vec fxx2(4);
arma::vec fxx3(4);
arma::vec fxx4(4);
arma::vec fxx5(4);
arma::vec fxx6(4);
arma::vec fxx7(4);
arma::vec fxx8(4);
arma::vec fxx9(4);
arma::vec fxx10(4);
arma::vec fxx11(4);
arma::vec fxx12(4);
arma::vec fxx13(4);
arma::vec fxx14(4);
arma::vec fxx15(4);
arma::vec fxx16(4);
arma::vec xx0(4);
arma::vec xx1(4);
arma::vec xx2(4);
arma::vec xx3(4);
arma::vec xx4(4);
arma::vec xx5(4);
arma::vec xx6(4);
arma::vec xx7(4);
arma::vec xx8(4);
arma::vec xx9(4);
arma::vec xx10(4);
arma::vec xx11(4);
arma::vec xx12(4);
arma::vec xx13(4);
arma::vec xx14(4);
arma::vec xx15(4);
arma::vec xx16(4);
xx0.fill(0);
xx0(0)=Bt;
arma::vec res1(N);
arma::vec res2(N);

double valA;
double valB;
double smtemp;
double d=tt;
for (int i = 1; i < N; i++)
{

fx0=f(x0,theta,d);
x1=x0+delt*(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
x10=x0+delt*(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
x11=x0+delt*(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
x12=x0+delt*(0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
x13=x0+delt*(0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15);
fx16=f(x16,theta,d+delt);

x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt;

valA=pcurr(x0,Bt,theta);

fxx0=f(xx0,theta,d);
xx1=xx0+delt*(0.1*fxx0);
fxx1=f(xx1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
xx2=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1);
fxx2=f(xx2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
xx3=xx0+delt*(0.202259190301118*fxx0+0.606777570903354*fxx2);
fxx3=f(xx3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt);
xx4=xx0+delt*(0.184024714708644*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3);
fxx4=f(xx4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
xx5=xx0+delt*(0.0879007340206681*fxx0+0.410459702520261*fxx3+0.482713753678866*fxx4);
fxx5=f(xx5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt);
xx6=xx0+delt*(0.085970050490246*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-0.0731856375070851*fxx5);
fxx6=f(xx6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
xx7=xx0+delt*(0.120930449125334*fxx0+0.260124675758296*fxx4+0.0325402621549091*fxx5-0.0595780211817361*fxx6);
fxx7=f(xx7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt);
xx8=xx0+delt*(0.110854379580391*fxx0-0.0605761488255006*fxx5+0.321763705601778*fxx6+0.510485725608063*fxx7);
fxx8=f(xx8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt);
xx9=xx0+delt*(0.112054414752879*fxx0-0.144942775902866*fxx5-0.333269719096257*fxx6+0.49926922955688*fxx7+0.509504608929686*fxx8);
fxx9=f(xx9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt);
xx10=xx0+delt*(0.113976783964186*fxx0-0.0768813364203357*fxx5+0.239527360324391*fxx6+0.397774662368095*fxx7+0.0107558956873607*fxx8-0.327769124164019*fxx9);
fxx10=f(xx10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt);
xx11=xx0+delt*(0.0798314528280196*fxx0-0.0520329686800603*fxx5-0.0576954146168549*fxx6+0.194781915712104*fxx7+0.145384923188325*fxx8-0.0782942710351671*fxx9-0.114503299361099*fxx10);
fxx11=f(xx11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt);
xx12=xx0+delt*(0.985115610164857*fxx0+0.330885963040722*fxx3+0.48966295730945*fxx4-1.37896486574844*fxx5-0.861164195027636*fxx6+5.78428813637537*fxx7+3.28807761985104*fxx8-2.38633905093136*fxx9-3.25479342483644*fxx10-2.16343541686423*fxx11);
fxx12=f(xx12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt);
xx13=xx0+delt*(0.895080295771633*fxx0+0.197966831227192*fxx2-0.0729547847313633*fxx3-0.851236239662008*fxx5+0.398320112318533*fxx6+3.63937263181036*fxx7+1.5482287703983*fxx8-2.12221714704054*fxx9-1.58350398545326*fxx10-1.71561608285936*fxx11-0.0244036405750127*fxx12);
fxx13=f(xx13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt);
xx14=xx0+delt*(-0.915176561375291*fxx0+1.45453440217827*fxx1+0*fxx2+0*fxx3-0.777333643644968*fxx4+0*fxx5-0.0910895662155176*fxx6+0.0910895662155176*fxx12+0.777333643644968*fxx13);
fxx14=f(xx14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt);
xx15=xx0+delt*(0.1*fxx0-0.157178665799771*fxx2+0.157178665799771*fxx14);
fxx15=f(xx15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt);
xx16=xx0+delt*(0.181781300700095*fxx0+0.675*fxx1+0.34275815984719*fxx2+0*fxx3+0.259111214548323*fxx4-0.358278966717952*fxx5-1.04594895940883*fxx6+0.930327845415627*fxx7+1.77950959431708*fxx8+0.1*fxx9-0.282547569539044*fxx10-0.159327350119973*fxx11-0.145515894647002*fxx12-0.259111214548323*fxx13-0.34275815984719*fxx14-0.675*fxx15);
fxx16=f(xx16,theta,d+delt);

xx0=xx0+(0.0333333333333333333333333333333333333333333333333333333333333*fxx0
+0.0250000000000000000000000000000000000000000000000000000000000*fxx1
+0.0333333333333333333333333333333333333333333333333333333333333*fxx2
+0.000000000000000000000000000000000000000000000000000000000000*fxx3
+0.0500000000000000000000000000000000000000000000000000000000000*fxx4
+0.000000000000000000000000000000000000000000000000000000000000*fxx5
+0.0400000000000000000000000000000000000000000000000000000000000*fxx6
+0.000000000000000000000000000000000000000000000000000000000000*fxx7
+0.189237478148923490158306404106012326238162346948625830327194*fxx8
+0.277429188517743176508360262560654340428504319718040836339472*fxx9
+0.277429188517743176508360262560654340428504319718040836339472*fxx10
+0.189237478148923490158306404106012326238162346948625830327194*fxx11
-0.0400000000000000000000000000000000000000000000000000000000000*fxx12
-0.0500000000000000000000000000000000000000000000000000000000000*fxx13
-0.0333333333333333333333333333333333333333333333333333333333333*fxx14
-0.0250000000000000000000000000000000000000000000000000000000000*fxx15
+0.0333333333333333333333333333333333333333333333333333333333333*fxx16)*delt;

d=d+delt;
valB=pcurr(xx0,Bt,theta);

res1(i)=valB;
if(i>1)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+res2(j-1)*res1(i-j);
}
res2(i) = (2.0*valA-2.0*smtemp*delt);
}
}
return(res2);
}

'

if(sum(theta[c(3,5,6)] ==0)==3)
{
cat(stre, " \r")
flush.console()
sourceCpp(code=txt1)
cat('                                     ','\r')
}
if(sum(theta[c(3,5,6)] ==0)!=3)
{
cat(stre, " \r")
flush.console()
sourceCpp(code=txt2)
cat('                                     ','\r')

}
N=round(t/delt)
res=solverFPT(Xs,B,c(0,theta),N,delt,0)
return(list(density=res, time = seq(0,t-delt,delt),prob.coverage = round(sum(res*delt),3)))
}
```

## Try the DiffusionRgqd package in your browser

Any scripts or data that you put into this service are public.

DiffusionRgqd documentation built on May 2, 2019, 3:26 a.m.