# # Case 9: # # In-line classical RK integrator # # h=0.1 - unstable; h=0.05 - stable # # ncase = 1: Classical fourth order Runge Kutta (RKC4) # # ncase = 2: Classical fourth order Runge Kutta with embedded # second order midpoint method and error estimate # # ODE routine bioreactor_9 (same as bioreactor_5, _6, _7, -8) # # No parms argument # # Derivative vector returned as numerical vector # rather than list # # ODE routine setwd("c:/R/bme_ode/chap1") source("bioreactor_9.R") # # Select classical fourth order Runge Kutta method # # ncase = 1: RKC4 # # ncase = 2: RKC4 with embedded second order midpoint method # and 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 # # rkc4 integration for(i1 in 2:nout){ # # nt rkc4 steps for(i2 in 1:nt){ if(ncase==1){ yb=y; tb=t rk1=bioreactor_9(tb,yb)*h y=yb+0.5*rk1; t=tb+0.5*h rk2=bioreactor_9(t,y)*h y=yb+0.5*rk2; t=tb+0.5*h rk3=bioreactor_9(t,y)*h y=yb+rk3; t=tb+h rk4=bioreactor_9(t,y)*h y=yb+(1/6)*(rk1+2*rk2+2*rk3+rk4) } if(ncase==2){ yb=y; tb=t rk1=bioreactor_9(tb,yb)*h y=yb+0.5*rk1; t=tb+0.5*h rk2=bioreactor_9(t,y)*h y2=yb+rk2 y=yb+0.5*rk2; t=tb+0.5*h rk3=bioreactor_9(t,y)*h y=yb+rk3; t=tb+h rk4=bioreactor_9(t,y)*h y=yb+(1/6)*(rk1+2*rk2+2*rk3+rk4) ee=y-y2 } } # # Solution after nt rkc4 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_9 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)