# # Rapid drug (anesthesia) delivery # # Access ODE integrator library("deSolve"); # # Access functions for numerical solutions setwd("c:/R/bme_pde/chap5"); source("pharma_2.R"); source("gs.R"); source("vf.R"); source("dss004.R"); source("dss012.R"); source("dss044.R"); # # Level of output # # ip = 1 - graphical (plotted) solutions # (u1(x,t), u2(x,t)) only # # ip = 2 - numerical and graphical solutions # ip=2; # # Grid in x nx=101;xl=0;xu=1; xg=seq(from=xl,to=xu,by=(xu-xl)/(nx-1)); # # Grid in t nout=51;t0=0;tf=10; tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1)); # # Parameters D1=3.0e-03;kb=1;u1e=1;ms=1; # # Display parameters cat(sprintf( "\n\n D1 = %8.3e kb = %6.3f u1e = %6.3f ms = %6.3f\n", D1,kb,u1e,ms)); Pe=vf(xg[51],0)*(xu-xl)/D1; cat(sprintf("\n Peclet number = %5.2f\n",Pe)); # # IC u0=rep(0,2*nx); ncall=0; # # ODE integration out=lsodes(y=u0,times=tout,func=pharma_2,parms=NULL) nrow(out) ncol(out) # # Arrays for plotting numerical solution u1_xplot=matrix(0,nrow=nx,ncol=nout); u2_xplot=matrix(0,nrow=nx,ncol=nout); for(it in 1:nout){ for(ix in 1:nx){ u1_xplot[ix,it]=out[it,ix+1]; u2_xplot[ix,it]=out[it,ix+1+nx]; } } # # Display numerical solution (for t = 0, 2,...,10) if(ip==2){ for(it in 1:nout){ if((it-1)*(it-11)*(it-21)*(it-31)*(it-41)*(it-51)==0){ cat(sprintf("\n\n t x u1(x,t) u2(x,t)")); for(ix in 1:nx){ cat(sprintf("\n%6.1f%7.3f%10.5f%10.5f", tout[it],xg[ix],u1_xplot[ix,it],u2_xplot[ix,it])); } } } } # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot u1, u2 numerical solutions # # vs x with t as a parameter, t = 0,0.2,...,10 par(mfrow=c(1,1)); matplot(x=xg,y=u1_xplot,type="l",xlab="x", ylab="u1(x,t), t=0,0.2,...,10",xlim=c(xl,xu),lty=1, main="u1(x,t); t=0,0.2,...,10;",lwd=2); par(mfrow=c(1,1)); matplot(x=xg,y=u1_xplot[,-1],type="l",xlab="x", ylab="u1(x,t), t=0.2,0.4,...,10",xlim=c(xl,xu),lty=1, main="u1(x,t); t=0.2,0.4,...,10;",lwd=2); par(mfrow=c(1,1)); matplot(x=xg,y=u2_xplot,type="l",xlab="x", ylab="u2(x,t), t=0,0.2,...,10",xlim=c(xl,xu),lty=1, main="u2(x,t); t=0,0.2,...,10;",lwd=2); # # vs t at x = 0.8, t = 0,0.2,...,10 par(mfrow=c(1,1)); u1_tplot=rep(0,nout); for(it in 1:nout){ u1_tplot[it]=u1_xplot[81,it]; } matplot(x=tout,y=u1_tplot,type="l",xlab="t", ylab="u1(x,t), x = 0.8",xlim=c(t0,tf),lty=1, main="u1(x,t); x = 0.8",lwd=2); par(mfrow=c(1,1)); u2_tplot=rep(0,nout); for(it in 1:nout){ u2_tplot[it]=u2_xplot[81,it]; } matplot(x=tout,y=u2_tplot,type="l",xlab="t", ylab="u2(x,t), x = 0.8",xlim=c(t0,tf),lty=1, main="u2(x,t); x = 0.8",lwd=2);