# # Corneal shape ODE/PDE # # Access deSolve library (with lsodes) library("deSolve") # # ODE/PDE routines setwd("c:/R/bme_ode/chap8"); source("corneal_1.R"); source("corneal_2.R"); source("dss006.R"); source("dss046.R"); # # Select case # # ncase = 1 - nonlinear ODE # # ncase = 2 - approximate nonlinear ODE # ncase=1; # # Parameters R=1; k=1; P=1; T=1; a=R^2*k/T; b=R*P/T; nx=21; xl=1; # # Write selected parameters cat(sprintf('\n\n xl = %5.2f a = %5.2f b = %5.2f\n\n',xl,a,b)); # # x grid, initial condition xg=seq(from=0,to=xl,by=xl/(nx-1)); u0=rep(0.25,nx); # # Output sequence in t tf=1;nout=6;ncall=0; tm=seq(from=0,to=tf,by=tf/(nout-1)); # # lsodes ODE integration parms=c(rtol=1e-8,atol=1e-8) out=lsodes(times=tm,y=u0,func=corneal_1,parms=parms); # # Numerical output cat(sprintf(" t u(0,t) u(0.25xl,t) u(0.50xl,t) u(0.75xl,t) u(xl,t)")); un=matrix(0,nrow=nout,ncol=nx); for(it in 1:nout){ if(it==1){un[1,]=u0; }else{ un[it,]=out[it,-1]; un[it,nx]=0;} cat(sprintf("%5.2f%13.5f%13.5f%13.5f%13.5f%13.5f\n", tm[it],un[it,1],un[it,6],un[it,11],un[it,16],un[it,21])); } # # Calls to corneal_1 cat(sprintf("\n ncall = %5d\n\n",ncall)) # # ux, ut for plotting ux_plot=matrix(0,nrow=nout,ncol=nx); ut_plot=matrix(0,nrow=nout,ncol=nx); for(it in 1:nout){ ux_plot[it,]=dss006(0,xl,nx,un[it,]); ut_plot[it,]=corneal_2(tm[it],un[it,],parms); } # # u(x,t) vs x matplot(x=xg,y=t(un),type="l",xlab="x",ylab="u(x,t)",xlim=c(0,xl), lty=1,main="u(x,t), t=0,0.2,...,1)",lwd = 2) # # ux(x,t) vs x matplot(x=xg,y=t(ux_plot),type="l",xlab="x",ylab="ux(x,t)",xlim=c(0,xl), lty=1,main="ux(x,t), t=0,0.2,...,1)",lwd = 2) # # ut(x,t) vs x matplot(x=xg,y=t(ut_plot[2:nout,]),type="l",xlab="x",ylab="ut(x,t)",xlim=c(0,xl), lty=1,main="ut(x,t), t=0.2,...,1)",lwd = 2)