# # Case 7: # # In-line modified Euler integrator # # h=0.1 - unstable; h=0.05 - stable # # ncase=1 - Runge Kutta format # # ncase=2 - ncase=1 with explicit error estimate # # ODE routine bioreactor_7 (same as bioreactor_5, _6) # # No parms argument # # Derivative vector returned as numerical vector # rather than list # # ODE routine setwd("c:/R/bme_ode/chap1") source("bioreactor_7.R") # # Select modified Euler format # # ncase = 1: RK format # # ncase = 2: RK format with explicit error estimate ncase=2; # # 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")) if(ncase==2){ ee=c(0,0,0,0,0,0,0) cat(sprintf( "\n e1 e2 e3 e4 e5 e6 e7"))} 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])) if(ncase==2){ cat(sprintf("\n %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n", ee[1],ee[2],ee[3],ee[4],ee[5],ee[6],ee[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 # # Modified Euler integration for(i1 in 2:nout){ # # nt modified Euler steps for(i2 in 1:nt){ if(ncase==1){ yb=y rk1=bioreactor_7(t,y)*h y=yb+rk1; t=t+h rk2=bioreactor_7(t,y)*h y=yb+(rk1+rk2)/2 } if(ncase==2){ yb=y rk1=bioreactor_7(t,y)*h y1=yb+rk1; t=t+h rk2=bioreactor_7(t,y1)*h y=yb+(rk1+rk2)/2 ee=y-y1 } } # # Solution after nt modified 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])) if(ncase==2){ cat(sprintf("\n %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n", ee[1],ee[2],ee[3],ee[4],ee[5],ee[6],ee[7]))} out[i1,-1]=y out[i1,1]=t } # # Calls to bioreactor_7 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)