# # Pattern formation model # # Murray, J.D. (2003), Mathematical Biology, II: Spatial Models # and Biomedical Applications, Springer, NY, p268 # # Access ODE integrator library("deSolve"); # # Access functions for numerical solutions setwd("c:/R/bme_pde/chap2/v_2pde"); source("p_form_1.R"); source("dss004.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; # # Initial condition (IC) # # ncase = 1 - spatially uniform # # ncase = 2 - Gaussian # # ncase = 3 - step # ncase=1; # # Grid (in x) nx=51;xl=0;xu=1; xg=seq(from=xl,to=xu,by=0.02); # # Parameters alpha=10;d1=0.5;w1=1;mu=0.1; cat(sprintf( "\n\n alpha = %5.2f d1 = %5.2f w1 = %5.2f mu = %5.2f\n", alpha,d1,w1,mu)); # # Independent variable for ODE integration nout=11; tout=seq(from=0,to=0.25,by=0.025); # # Initial condition u0=rep(0,2*nx);u10=rep(0,nx);u20=rep(0,nx); for(i in 1:nx){ # # Solution remains spatially uniform if(ncase==1){ u10[i]=1;u20[i]=0;} # # Initial Gaussian in u1 if(ncase==2){ u10[i]=exp(-5*(xg[i]-0.5)^2);u20[i]=0;} # # Initial band in u1 if(ncase==3){ if(i<=11){ u10[i]=1;u20[i]=0; }else{ u10[i]=0;u20[i]=0; } } u0[i] =u10[i]; u0[i+nx]=u20[i]; } t=0; ncall=0; # # ODE integration out=lsodes(y=u0,times=tout,func=p_form_1,parms=NULL) nrow(out) ncol(out) # # Arrays for plotting numerical solution u1_plot=matrix(0,nrow=nx,ncol=nout); u2_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]; } } # # Display numerical solution if(ip==2){ for(it in 1:nout){ cat(sprintf("\n t x u1(x,t) u2(x,t)\n")); for(ix in 1:nx){ cat(sprintf("%7.3f%8.2f%12.5f%12.5f\n", tout[it],xg[ix],u1_plot[ix,it],u2_plot[ix,it])); } } } # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot u1 numerical par(mfrow=c(1,1)); matplot(x=xg,y=u1_plot,type="l",xlab="x", ylab="u1(x,t), t=0,0.025,...,0.25",xlim=c(xl,xu),lty=1, main="u1(x,t); t=0,0.025,...,0.25;",lwd=2); # # Plot u2 numerical par(mfrow=c(1,1)); matplot(x=xg,y=u2_plot,type="l",xlab="x", ylab="u2(x,t), t=0,0.025,...,0.25",xlim=c(xl,xu),lty=1, main="u2(x,t); t=0,0.025,...,0.25;",lwd=2);