# # Fitzhugh-Nagumo (FHN) equation # # Access ODE integrator library("deSolve"); # # Access functions for numerical solution setwd("c:/R/bme_pde/chap4"); source("fhn_1.R"); source("dss004.R"); # # Level of output # # ip = 1 - graphical (plotted) solutions # (u(x,t)) only # # ip = 2 - numerical and graphical solutions # ip=2; # # Alternative boundary conditions (BCs) # # ncase = 1 - analytical Dirichlet BCs # (not used) # # ncase = 2 - constant Dirichlet BCs # # ncase = 3 - analytical Neumann BCs # (not used) # # ncase = 4 - homogeneous Neumann BCs # # ncase = 5 - no BCs # ncase=2; # # Parameters a=1;D=1; cat(sprintf("\n\n a = %5.2f D = %5.2f\n",a,D)); # # Grid (in x) nx=101;xl=-60;xu=20; xg=seq(from=xl,to=xu,by=(xu-xl)/(nx-1)); # # Independent variable for ODE integration nout=7; tout=seq(from=0,to=60,by=10); # # Initial condition (unit step) u0=rep(0,nx);t0=0; for(i in 1:101){ if(i< 76){u0[i]=1; } if(i==76){u0[i]=0.5;} if(i> 76){u0[i]=0; } } ncall=0; # # ODE integration out=lsodes(y=u0,times=tout,func=fhn_1,parms=NULL) nrow(out) ncol(out) # # Arrays for plotting numerical solution u_plot=matrix(0,nrow=nx,ncol=nout); for(it in 1:nout){ for(ix in 1:nx){ u_plot[ix,it]=out[it,ix+1]; } } # # Display numerical solution if(ip==2){ for(it in 1:nout){ cat(sprintf("\n t x u(x,t)\n")); for(ix in 1:nx){ cat(sprintf("%5.1f%8.2f%10.5f\n",tout[it],xg[ix],u_plot[ix,it])); } } } # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot u numerical, analytical solution par(mfrow=c(1,1)); matplot(x=xg,y=u_plot,type="l",xlab="x",ylab="u1(x,t),t=0,10,...,60", xlim=c(xl,xu),lty=1,main="u(x,t);t=0,10,...,60;",lwd=2);