# # Case 5: # # In-line Euler integrator # # h=0.1, unstable # # h=0.05; stable # # ODE routine bioreactor_5 # # No parms argument # # Derivative vector returned as numerical vector # rather than list # # ODE routine setwd("c:/R/bme_ode/chap1") source("bioreactor_5.R") # # Parameter values for BP10001 k1=8.87e-03; k2=13.18; k3=0.129; k4=0.497; k5=0.027; k6=0.545e-3; km2=87.7; km3=99.9; # # Initial condition n=7;nout=51;t=0;ncall=0; y=c(0.10724,0,0,0,0,0,0) cat(sprintf( "\n t y1 y2 y3 y4 y5 y6 y7")) cat(sprintf("\n %8.0f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f", t,y[1],y[2],y[3],y[4],y[5],y[6],y[7])) # # Arrays for output out=matrix(0,nrow=nout,ncol=(n+1)) out[1,-1]=y out[1,1]=t # # Parameters for t integration # nt=400;h=0.10 # unstable nt=800;h=0.05 # stable # # Euler integration for(i1 in 2:nout){ # # nt Euler steps for(i2 in 1:nt){ yt=bioreactor_5(t,y); y=y+yt*h; t=t+h; } # # Solution after nt Euler steps cat(sprintf("\n %8.0f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f", t,y[1],y[2],y[3],y[4],y[5],y[6],y[7])) out[i1,-1]=y out[i1,1]=t } # # Calls to bioreactor_5 cat(sprintf("\n ncall = %5d\n\n",ncall)) # # Single plot par(mfrow=c(1,1)) # # y1 plot(out[,1],out[,2],type="l",xlab="t",ylab="y1(t),...,y7(t)", xlim=c(0,2000),ylim=c(0,0.14),lty=1, main="y1(t),...,y7(t) vs t", lwd=2) # # y2 lines(out[,1],out[,3],type="l",lty=2,lwd=2) # # y3 lines(out[,1],out[,4],type="l",lty=3,lwd=2) # # y4 lines(out[,1],out[,5],type="l",lty=4,lwd=2) # # y5 lines(out[,1],out[,6],type="l",lty=5,lwd=2) # # y6 lines(out[,1],out[,7],type="l",lty=6,lwd=2) # # y7 lines(out[,1],out[,8],type="l",lty=7,lwd=2)