pde_terms=function(t,u){ # # Function pde_terms provides supplmental calculations for the three-pde # model based on the u1, u2, u3 vectors # # One vector to three vectors u1=rep(0,nx);u2=rep(0,nx);u3=rep(0,nx); for(i in 1:nx){ u1[i]=u[i]; u2[i]=u[i+nx]; u3[i]=u[i+2*nx]; } # # u1x, u2x u1x=dss004(xl,xu,nx,u1); u2x=dss004(xl,xu,nx,u2); u3x=dss004(xl,xu,nx,u3); # # Boundary conditions u1x[1]=0;u1x[nx]=0; u2x[1]=0;u2x[nx]=0; u3x[1]=0;u3x[nx]=0; nl=2;nu=2; # # u1xx, u2xx u1xx=dss044(xl,xu,nx,u1,u1x,nl,nu); u2xx=dss044(xl,xu,nx,u2,u2x,nl,nu); u3xx=dss044(xl,xu,nx,u3,u3x,nl,nu); # # RHS terms term1=rep(0,nx);term2=rep(0,nx);term3=rep(0,nx); for(i in 1:nx){ den=1/(k[2]+u2[i])^2; term1[i]=k[1]*u1[i]*den*u2xx[i]; term2[i]=k[1]*den*u1x[i]*u2x[i]; term3[i]=-2*k[1]*u1[i]*den/(k[2]+u2[i])*u2x[i]^2; } # # PDEs u1t=rep(0,nx);u2t=rep(0,nx);u3t=rep(0,nx); for(i in 1:nx){ u1t[i]=D1*u1xx[i]-(term1[i]+term2[i]+term3[i])+ k[3]*u1[i]*(k[4]*u3[i]^2/(k[9]+u3[i]^2)-u1[i]); u2t[i]=D2*u2xx[i]+k[5]*u3[i]*(u1[i]^2/(k[6]+u1[i]^2)- k[7]*u1[i]*u2[i]); u3t[i]=D3*u3xx[i]-k[8]*u1[i]*(u3[i]^2/(k[9]+u3[i]^2)); } # # PDE terms chemo1=rep(0,nx);chemo2=rep(0,nx);chemo3=rep(0,nx); chemo4=rep(0,nx);chemo5=rep(0,nx);chemo6=rep(0,nx); chemo7=rep(0,nx);chemo8=rep(0,nx);chemo9=rep(0,nx); chemo10=rep(0,nx); for(i in 1:nx){ chemo1[i]=D1*u1xx[i]; chemo2[i]=-(term1[i]+term2[i]+term3[i]); chemo3[i]=k[3]*u1[i]*(k[4]*u3[i]^2/(k[9]+u3[i]^2)-u1[i]); chemo4[i]=D2*u2xx[i]; chemo5[i]=k[5]*u3[i]*(u1[i]^2/(k[6]+u1[i]^2)-k[7]*u1[i]*u2[i]); chemo6[i]=D3*u3xx[i]; chemo7[i]=-k[8]*u1[i]*(u3[i]^2/(k[9]+u3[i]^2)); chemo8[i]=u1t[i]; chemo9[i]=u2t[i]; chemo10[i]=u3t[i]; } # # Return arrays of PDE terms chemo1 <<- chemo1;chemo2 <<- chemo2;chemo3 <<- chemo3; chemo4 <<- chemo4;chemo5 <<- chemo5;chemo6 <<- chemo6; chemo7 <<- chemo7;chemo8 <<- chemo8;chemo9 <<- chemo9; chemo10<<-chemo10; }