# # Delete previous workspaces rm(list=ls(all=TRUE)) # # Access ODE integrator library("deSolve"); # # Access ODE file setwd("g:/matlab_R/ode"); source("ode_1.R"); # # Matrices for plotted output nout=21; Vplot=matrix(0,nrow=nout,ncol=3); Vaplot=matrix(0,nrow=nout,ncol=3); diff=matrix(0,nrow=nout,ncol=3); # # Initial condition lambda=rep(0,3); for(icase in 1:3){ if(icase==1){lambda[1]=1 ;} if(icase==2){lambda[2]=0.75;} if(icase==3){lambda[3]=1.25;} V0=1;alpha=1;ncall=0; # # t scale tout=seq(from=0,to=10,by=0.5); # # ODE integration out=ode(func=ode_1,times=tout,y=V0,parms); # # Display heading V=rep(0,nout);Va=rep(0,nout); V=out[,-1]; cat(sprintf("\n\n icase = %2d\n",icase)); # # Display calls to ODE routine cat(sprintf("\n ncall = %3d\n",ncall)); # # Display numerical output cat(sprintf("\n t V(t)(num) V(t)(anal) diff\n")); for(it in 1:nout){ Va[it]=V0*exp(lambda[icase]/alpha*(1-exp(-alpha*tout[it]))); Vplot[it,icase]= V[it]; Vaplot[it,icase]=Va[it]; diff[it,icase]=(Vplot[it,icase]-Vaplot[it,icase])/ Vplot[it,icase]; cat(sprintf("%10.2f%12.4f%12.4f%12.6f\n", tout[it],Vplot[it,icase],Vaplot[it,icase],diff[it,icase])); } } # # Plot numerical, analytical solutions par(mfrow=c(1,1)) plot(tout,Vplot[,1],xlab="t",ylab="V(t)", xlim=c(0,10),ylim=c(1,5), main="V(t), solid - num, points - anal",type="l",lwd=2); points(tout,Vaplot[,1],pch="1"); lines(tout, Vplot[,2],type="l",lwd=2); points(tout,Vaplot[,2],pch="2"); lines(tout, Vplot[,3],type="l",lwd=2); points(tout,Vaplot[,3],pch="3"); # # Plot difference of solutions par(mfrow=c(1,1)) matplot(tout,diff,type="l",xlab="t",ylab="diff", col="black",lwd=2,lty=1,main="diff = num - anal (rel)");