# # Drug tracking # # Access ODE integrator library("deSolve"); # # Access functions for numerical solutions setwd("c:/R/bme_pde/chap7"); source("drug_1.R"); source("drug_2.R"); source("f_u.R"); source("g_v.R"); source("dss004.R"); source("dss044.R"); # # Format of output # # ip = 1 - graphical (plotted) solutions vs x with # t as a parameter # # ip = 2 - graphical solutions vs t at specific x # ip=1; # # Select ODE routine # # ncase = 1 - stagewise differentiation # # ncase = 2 - direct second derivative # ncase=1; # # Grid in x nx=26;xl=-0.5;xu=0.5; xg=seq(from=xl,to=xu,by=(xu-xl)/(nx-1)); # # Grid in t if(ip==1){nout= 6;t0=0;tf=2;} if(ip==2){nout=41;t0=0;tf=2;} tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1)); # # Parameters alpha=0.2; beta=1; gamma=1; D=0.6; E=0.2; kr=1; ub=1; vb=1; ua=0; # # Display selected parameters cat(sprintf( "\n\n D = %6.3f E = %6.3f\n",D,E)); # # IC u0=rep(0,3*nx); for(ix in 1:nx){ u0[ix] =0.75; u0[ix+nx] =0.25; u0[ix+2*nx]=0; } ncall=0; # # ODE integration if(ncase==1){ out=ode(y=u0,times=tout,func=drug_1,parms=NULL);} if(ncase==2){ out=ode(y=u0,times=tout,func=drug_2,parms=NULL);} nrow(out) ncol(out) # # Arrays for plotting numerical solutions u_xplot=matrix(0,nrow=nx,ncol=nout); v_xplot=matrix(0,nrow=nx,ncol=nout); s_xplot=matrix(0,nrow=nx,ncol=nout); for(it in 1:nout){ for(ix in 1:nx){ u_xplot[ix,it]=out[it,ix+1]; v_xplot[ix,it]=out[it,ix+1+nx]; s_xplot[ix,it]=out[it,ix+1+2*nx]; } } # # Display numerical solutions (for t = 0, 2) if(ip==1){ for(it in 1:nout){ if((it-1)*(it-6)==0){ cat(sprintf( "\n\n t x u(x,t) v(x,t) s(x,t)")); for(ix in 1:nx){ cat(sprintf("\n %6.1f%7.2f%12.5f%12.5f%12.5f", tout[it],xg[ix],u_xplot[ix,it],v_xplot[ix,it],s_xplot[ix,it])); } } } } if(ip==2){ for(it in 1:nout){ if((it-1)*(it-41)==0){ cat(sprintf( "\n\n t x u(x,t) v(x,t) s(x,t)")); for(ix in 1:nx){ cat(sprintf("\n %6.1f%7.2f%12.5f%12.5f%12.5f", tout[it],xg[ix],u_xplot[ix,it],v_xplot[ix,it],s_xplot[ix,it])); } } } } # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot u,v,s numerical solutions # # vs x with t as a parameter, t = 0,0.4,...,2 if(ip==1){ # # t=0 par(mfrow=c(1,1)); plot(xg,u_xplot[,1],type="b",lty=1,pch="1",xlab="x", ylab="u(x,t),t=0,0.4,...,2",xlim=c(xl,xu),ylim=c(0,0.8), main="u(x,t); t=0,0.4,...,2;",lwd=2) # # t=0.4,0.8,1.2,1.6,2 lines(xg,u_xplot[,2],type="b",lty=1,pch="2",lwd=2) lines(xg,u_xplot[,3],type="b",lty=1,pch="3",lwd=2) lines(xg,u_xplot[,4],type="b",lty=1,pch="4",lwd=2) lines(xg,u_xplot[,5],type="b",lty=1,pch="5",lwd=2) lines(xg,u_xplot[,6],type="b",lty=1,pch="6",lwd=2) # # t=0 par(mfrow=c(1,1)); plot(xg,v_xplot[,1],type="b",lty=1,pch="1",xlab="x", ylab="v(x,t),t=0,0.4,...,2",xlim=c(xl,xu),ylim=c(0.1,0.3), main="v(x,t); t=0,0.4,...,2;",lwd=2) # # t=0.4,0.8,1.2,1.6,2 lines(xg,v_xplot[,2],type="b",lty=1,pch="2",lwd=2) lines(xg,v_xplot[,3],type="b",lty=1,pch="3",lwd=2) lines(xg,v_xplot[,4],type="b",lty=1,pch="4",lwd=2) lines(xg,v_xplot[,5],type="b",lty=1,pch="5",lwd=2) lines(xg,v_xplot[,6],type="b",lty=1,pch="6",lwd=2) # # t=0 par(mfrow=c(1,1)); plot(xg,s_xplot[,1],type="b",lty=1,pch="1",xlab="x", ylab="s(x,t),t=0,0.4,...,2",xlim=c(xl,xu),ylim=c(-0.4,0.0), main="s(x,t); t=0,0.4,...,2;",lwd=2) # # t=0.4,0.8,1.2,1.6,2 lines(xg,s_xplot[,2],type="b",lty=1,pch="2",lwd=2) lines(xg,s_xplot[,3],type="b",lty=1,pch="3",lwd=2) lines(xg,s_xplot[,4],type="b",lty=1,pch="4",lwd=2) lines(xg,s_xplot[,5],type="b",lty=1,pch="5",lwd=2) lines(xg,s_xplot[,6],type="b",lty=1,pch="6",lwd=2) } # # vs t at x = 0, t = 0,0.05,...,2 if(ip==2){ u_tplot=rep(0,nout);v_tplot=rep(0,nout); s_tplot=rep(0,nout); for(it in 1:nout){ u_tplot[it]=u_xplot[13,it]; v_tplot[it]=v_xplot[13,it]; s_tplot[it]=s_xplot[13,it]; } par(mfrow=c(1,1)); matplot(x=tout,y=u_tplot,type="l",xlab="t", ylab="u(x,t), x = 0",xlim=c(t0,tf),lty=1, main="u(x,t); x = 0",lwd=2); par(mfrow=c(1,1)); matplot(x=tout,y=v_tplot,type="l",xlab="t", ylab="v(x,t), x = 0",xlim=c(t0,tf),lty=1, main="v(x,t); x = 0",lwd=2); par(mfrow=c(1,1)); matplot(x=tout,y=s_tplot,type="l",xlab="t", ylab="s(x,t), x = 0",xlim=c(t0,tf),lty=1, main="s(x,t); x = 0",lwd=2); }