MOL.passage=function(Xs,t,limits,N,delt,mu,sig,desc=1)
{
xx1=seq(limits[1],limits[2],length=N)
dx1=diff(xx1)[1]
dxx1=dx1^2
if((missing(mu)&&missing(sig)))
{
namess=c('mu','sig')
namess2=c('muu','sigg')
txt =rep('+rep(0,length(X))',2)
func.list=rep(0,length(namess))
obs=objects(pos=1)
muu =function(X,t){}
sigg =function(X,t){}
for(i in 1:length(namess))
{
if(sum(obs==namess[i]))
{
txt[i] = paste0(body(namess[i])[2],txt[i])
}
}
body(muu) = parse(text=txt[1])
body(sigg) = parse(text=txt[2])
}else
{
muu =function(X,t){}
sigg =function(X,t){}
body(muu) = parse(text=paste0(mu,'+rep(0,length(X))'))
body(sigg)= parse(text=paste0(sig,'+rep(0,length(X))'))
}
if(desc==1)
{
f=function(U,tme)
{
SU1= sigg(xx1)^2
D1 = muu(xx1)[-c(1,N)]*(U[-c(1,2)]-U[-c(N-1,N)])/dx1/2
D2 = 1/2*SU1[-c(1,N)]*((U[-c(1,2)]-2*U[-c(1,N)]+U[-c(N-1,N)]))/dxx1
MMM1=(D1+D2)
return(c(0,MMM1,0))
}
M=rep(1,N)
M[c(1,N)]=0
}
if(desc==2)
{
f=function(U,tme)
{
SU1= sigg(xx1)^2
D1 = muu(xx1)[-N]*(U[-1]-U[-N])/dx1/2
D2 = 1/2*SU1[-c(1,N)]*((U[-c(1,2)]-2*U[-c(1,N)]+U[-c(N-1,N)]))/dxx1
MMM1=(D1+c(0,D2))
return(c(MMM1,0))
}
M=rep(1,N)
M[c(N)]=0
}
if(desc==3)
{
f=function(U,tme)
{
SU1= sigg(xx1)^2
D1 = muu(xx1)[-1]*(U[-1]-U[-N])/dx1/2
D2 = 1/2*SU1[-c(1,N)]*((U[-c(1,2)]-2*U[-c(1,N)]+U[-c(N-1,N)]))/dxx1
MMM1=(D1+c(D2,0))
return(c(0,MMM1))
}
M=rep(1,N)
M[c(1)]=0
}
N.mesh=round((t-0)/delt)+1
MM = array(0,dim=c(N,N.mesh))
MM[,1]=M
ttt=0
t.alpha=
c(0.000000000000000000000000000000000000000000000000000000000000,
0.100000000000000000000000000000000000000000000000000000000000,
0.539357840802981787532485197881302436857273449701009015505500,
0.809036761204472681298727796821953655285910174551513523258250,
0.309036761204472681298727796821953655285910174551513523258250,
0.981074190219795268254879548310562080489056746118724882027805,
0.833333333333333333333333333333333333333333333333333333333333,
0.354017365856802376329264185948796742115824053807373968324184,
0.882527661964732346425501486979669075182867844268052119663791,
0.642615758240322548157075497020439535959501736363212695909875,
0.357384241759677451842924502979560464040498263636787304090125,
0.117472338035267653574498513020330924817132155731947880336209,
0.833333333333333333333333333333333333333333333333333333333333,
0.309036761204472681298727796821953655285910174551513523258250,
0.539357840802981787532485197881302436857273449701009015505500,
0.100000000000000000000000000000000000000000000000000000000000,
1.00000000000000000000000000000000000000000000000000000000000)
for(i in 2:N.mesh)
{
x0=M
fx0=f(x0,ttt)
x1=x0+delt*(0.1*fx0)
fx1=f(x1,ttt+t.alpha[2]*delt)
x2=x0+delt*(-0.915176561375291*fx0 +1.45453440217827*fx1)
fx2=f(x2,ttt+t.alpha[3]*delt)
x3=x0+delt*( 0.202259190301118*fx0 +0.606777570903354*fx2)
fx3=f(x3,ttt+t.alpha[4]*delt)
x4=x0+delt*( 0.184024714708644*fx0 +0.197966831227192*fx2-0.0729547847313633*fx3)
fx4=f(x4,ttt+t.alpha[5]*delt)
x5=x0+delt*( 0.0879007340206681*fx0 +0.410459702520261*fx3+0.482713753678866*fx4)
fx5=f(x5,ttt+t.alpha[6]*delt)
x6=x0+delt*(0.085970050490246*fx0 +0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5)
fx6=f(x6,ttt+t.alpha[7]*delt)
x7=x0+delt*(0.120930449125334*fx0 +0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6)
fx7=f(x7,ttt+t.alpha[8]*delt)
x8=x0+delt*(0.110854379580391*fx0 -0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7)
fx8=f(x8,ttt+t.alpha[9]*delt)
x9=x0+delt*(0.112054414752879*fx0 -0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8)
fx9=f(x9,ttt+t.alpha[10]*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,ttt+t.alpha[11]*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,ttt+t.alpha[12]*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,ttt+t.alpha[13]*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,ttt+t.alpha[14]*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,ttt+t.alpha[15]*delt)
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14)
fx15=f(x15,ttt+t.alpha[16]*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,ttt+t.alpha[17]*delt)
M=M+(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
ttt=ttt+delt
MM[,i]=M
if(any(is.na(M))){stop(paste("NAs produced at time",ttt,". Perhaps increase NN and decrease delt or change desc."))}
}
if(sum(abs(xx1-Xs[1])==0)==1)
{
wh1=which(abs(xx1-Xs[1])==0)
dens=c(0,-diff(MM[wh1,])/delt)
}
if(sum(abs(xx1-Xs[1])==0)==0)
{
tst1=(Xs-xx1)
wh1=min(which(tst1<0))
wh2=max(which(tst1>=0))
ymin=MM[wh1,]
ymax=MM[wh2,]
y=ymin+(ymax-ymin)*(Xs-xx1[wh2])/(xx1[wh1]-xx1[wh2])
dens=c(0,-diff(y)/delt)
}
ret = (list(surface=MM,density=dens,Xt=xx1,time=seq(0,t,delt)))
class(ret) = 'MOL.passage'
return(ret)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.