p_form_2=function(t,u,parms){ # # Function p_form_2 computes the t derivative vector of 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, u3x 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, u3xx 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)); } # # Three vectors to one vector ut=rep(0,3*nx); for(i in 1:nx){ ut[i] =u1t[i]; ut[i+nx] =u2t[i]; ut[i+2*nx]=u3t[i]; } # # Increment calls to p_form_2 ncall <<- ncall+1; # # Return derivative vector return(list(c(ut))); }