# # Two PDE chemotaxis model # # Murray, J.D. (2003), Mathematical Biology, II: Spatial Models # and Biomedical Applications, Springer, NY, p68, prob. 5 # # Access ODE integrator library("deSolve"); # # Access functions for numerical, analytical solutions setwd("c:/R/bme_Rroutines/temp"); source("chemo_1.R"); source("u1_anal.R"); source("u2_anal.R"); source("dss004.R"); # # Level of output # # ip = 1 - graphical (plotted) solutions # (u1(x,t), u2(x,t)) only # # ip = 2 - numerical and graphical solutions # # ip = 3 - ip = 1 plus graphical output of # PDE RHS terms # ip=3; # # Grid (in x) nx=101;xl=-10;xu=15 xg=seq(from=xl,to=xu,by=0.25); # # Parameters k=1;D=1;c=1; cat(sprintf("\n\n k = %5.2f D = %5.2f c = %5.2f\n",k,D,c)); # # Independent variable for ODE integration nout=6; tout=seq(from=0,to=5,by=1); # # Initial condition (from analytical solutions,t=0) u0=rep(0,2*nx); for(i in 1:nx){ u0[i] =u1_anal(xg[i],tout[1],k,D,c); u0[i+nx]=u2_anal(xg[i],tout[1],k,D,c); } ncall=0; # # ODE integration out=lsodes(y=u0,times=tout,func=chemo_1,parms=NULL) nrow(out) ncol(out) # # Arrays for plotting numerical, analytical solutions u1_plot=matrix(0,nrow=nx,ncol=nout); u2_plot=matrix(0,nrow=nx,ncol=nout); u1a_plot=matrix(0,nrow=nx,ncol=nout); u2a_plot=matrix(0,nrow=nx,ncol=nout); for(it in 1:nout){ for(ix in 1:nx){ u1_plot[ix,it]=out[it,ix+1]; u2_plot[ix,it]=out[it,ix+1+nx]; u1a_plot[ix,it]=u1_anal(xg[ix],tout[it],k,D,c); u2a_plot[ix,it]=u2_anal(xg[ix],tout[it],k,D,c); } } # # Display numerical solution if(ip==2){ for(it in 1:nout){ cat(sprintf("\n t x u1(x,t) u1_ex(x,t) u1_err(x,t)")); cat(sprintf("\n u2(x,t) u2_ex(x,t) u2_err(x,t)\n")); for(ix in 1:nx){ cat(sprintf("%5.1f%8.2f%10.5f%12.5f%13.6f\n",tout[it],xg[ix], u1_plot[ix,it],u1a_plot[ix,it],u1_plot[ix,it]-u1a_plot[ix,it])); cat(sprintf(" %10.5f%12.5f%13.6f\n", u2_plot[ix,it],u2a_plot[ix,it],u2_plot[ix,it]-u2a_plot[ix,it])); } } } # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot u1 numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u1_plot,type="l",xlab="x",ylab="u1(x,t),t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u1(x,t); solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u1a_plot,xlim=c(xl,xu),col="black",lwd=2) # # Plot u2 numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u2_plot,type="l",xlab="x",ylab="u2(x,t),t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u2(x,t), solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u2a_plot,xlim=c(xl,xu),col="black",lwd=2) # # Plotting of PDE RHS terms if(ip==3){ # # 1D arrays of various PDE RHS terms (denoted 1d) u1_1d =rep(0,nx); u2_1d =rep(0,nx);u1u2_1d=rep(0,nx); u1x_1d =rep(0,nx);u2x_1d =rep(0,nx); u1u2x_1d=rep(0,nx);u2xx_1d=rep(0,nx); # # 2D arrays for plotting (denoted 2d) u1x_2d =matrix(0,nrow=nx,ncol=nout); u2x_2d =matrix(0,nrow=nx,ncol=nout); u1u2_2d =matrix(0,nrow=nx,ncol=nout); u2xx_2d =matrix(0,nrow=nx,ncol=nout); u1u2x_2d=matrix(0,nrow=nx,ncol=nout); u1t_2d =matrix(0,nrow=nx,ncol=nout); u2t_2d =matrix(0,nrow=nx,ncol=nout); u1x_2d_anal=matrix(0,nrow=nx,ncol=nout); u2x_2d_anal=matrix(0,nrow=nx,ncol=nout); u1t_2d_anal=matrix(0,nrow=nx,ncol=nout); u2t_2d_anal=matrix(0,nrow=nx,ncol=nout); # # PDE RHS terms for(it in 1:nout){ for(ix in 1:nx){ u1_1d[ix]=u1_plot[ix,it]; u2_1d[ix]=u2_plot[ix,it]; } u1x_1d=dss004(xl,xu,nx,u1_1d); u2x_1d=dss004(xl,xu,nx,u2_1d); u1x_1d[1]=0;u1x_1d[nx]=0; u2x_1d[1]=0;u2x_1d[nx]=0; for(ix in 1:nx){ u1u2_1d[ix]=u2_1d[ix]/u1_1d[ix]*u1x_1d[ix]; } u1u2x_1d=dss004(xl,xu,nx,u1u2_1d); u2xx_1d=dss004(xl,xu,nx, u2x_1d); # # 2D arrays for plotting for(ix in 1:nx){ # # Derivatives of solutions in x u1x_2d[ix,it]=u1x_1d[ix]; u2x_2d[ix,it]=u2x_1d[ix]; u1u2_2d[ix,it]=u1u2_1d[ix]; u2xx_2d[ix,it]=u2xx_1d[ix]; u1u2x_2d[ix,it]=u1u2x_1d[ix]; # # Derivatives of solutions in t u1t_2d[ix,it]=-k*u2_1d[ix]; u2t_2d[ix,it]=D*(u2xx_1d[ix]-2*u1u2x_1d[ix]); # # Analytical derivatives of solutions in t expz=exp(-c*(xg[ix]-c*tout[it])/D); u1x_2d_anal[ix,it] =(1/c)*(1/(1+expz)^2)*expz*(c^2/D); u2x_2d_anal[ix,it]=-(1/c)* (c^4/(k*D^2))*expz*(1-expz)/(1+expz)^3; u1t_2d_anal[ix,it]=-c*u1x_2d_anal[ix,it]; u2t_2d_anal[ix,it]=-c*u2x_2d_anal[ix,it]; } # # Next t } # # Plot Du2_{xx} par(mfrow=c(1,1)); matplot(x=xg,y=D*u2xx_2d,type="l",xlab="x",ylab="Du2_{xx},t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="Du2_{xx}; t=0,1,2,3,4,5;",lwd=2); # # Plot -2D((u2/u1)u1_x)_x par(mfrow=c(1,1)); matplot(x=xg,y=-2*D*u1u2x_2d,type="l",xlab="x",ylab="-2D((u2/u1)u1_x)_x, t=0,1,2,3,4,5",xlim=c(xl,xu),lty=1 ,main="-2D((u2/u1)u1_x)_x, t=0,1,2,3,4,5;",lwd=2); # # Plot u1 derivative in x numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u1x_2d,type="l",xlab="x",ylab="u1(x,t)_x,t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u1(x,t)_x; solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u1x_2d_anal,xlim=c(xl,xu),col="black",lwd=2) # # Plot u2 derivative in x numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u2x_2d,type="l",xlab="x",ylab="u2(x,t)_t,t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u2(x,t)_x, solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u2x_2d_anal,xlim=c(xl,xu),col="black",lwd=2) # # Plot u1 derivative in t numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u1t_2d,type="l",xlab="x",ylab="u1(x,t)_t,t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u1(x,t)_t; solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u1t_2d_anal,xlim=c(xl,xu),col="black",lwd=2) # # Plot u2 derivative in t numerical, analytical par(mfrow=c(1,1)); matplot(x=xg,y=u2t_2d,type="l",xlab="x",ylab="u2(x,t)_t,t=0,1,2,3,4,5", xlim=c(xl,xu),lty=1,main="u2(x,t)_t, solid - num, points - anal; t=0,1,2,3,4,5;",lwd=2); matpoints(x=xg,y=u2t_2d_anal,xlim=c(xl,xu),col="black",lwd=2) # # End ip = 3 }