2
$\begingroup$

I'm trying to do something very similar to this question, except with a much larger system of PDEs from this paper. The PDEs are equations 1-18 in that paper. Equations 9, 10 and 15-18 have laplacian terms in them. From the paper, I know that a rectangular domain is used and that the variables that have laplacian terms in their equations satisfy periodic boundary conditions. Therefore, I wrote the following code:

 (* Parameters *) DAO = 4.32*10^-2; DH = 8.11*10^-2; DTAlpha = 6.55*10^-2; DTBeta = 6.55*10^-2; Di10 = 6.04*10^-2; DP = 1.2*10^-1; lambdaBetaI = 9.51*10^-6; lambdaN = 8*10^-9; lambdaA = 8*10^-10; lambdaTau0 = 8.1*10^-11; lambdaTau = 1.35*10^-11; lambdaF = 1.662*10^-3; lambdaATAlpha = 1.54; lambdaAABeta0 = 1.793; lambdaAO = 5*10^-2; lambdaH = 3*10^-5; lambdaMF = 2*10^-2; lambdaMA = 2.3*10^-3; lambdaM1TBeta = 6*10^-3; lambdaM1HatTBeta = 6*10^-4; lambdaTBetaM2 = 1.5 * 10^-2; lambdaTBetaM2Hat = 1.5 * 10^-2; lambdaTAlphaM1 = 3 * 10^-2; lambdaTAlphaM1Hat = 3 * 10^-2; lambdai10M2 = 6.67 * 10^-3; lambdai10M2Hat = 6.67 * 10^-3; lambdaPA = 6.6 * 10^-8; lambdaPM2 = 1.32 * 10^-7; myTheta = 0.9; alpha1 = 5; myBeta = 10; myGamma = 1; dABetaI = 9.51; dABeta0M = 9.51; dABeta0M = 2 * 10^-3; dABeta0MHat = 10^-2; dTau = 0.277; dFI = 2.77 * 10^-3; dFO = 2.77 * 10^-4; dN = 1.9 * 10^-4; dNF = 3.4 * 10^-4; dNT = 1.7 * 10^-4; dNdM = 0.06; dNdMhat = 0.02; dA = 1.2 * 10^-3; dM1 = 0.015; dM2 = 0.015; dM1Hat = 0.015; dM2Hat = 0.015; dAO = 0.951; dH = 58.71; dTAlpha = 55.45; dTBeta = 3.33 * 10^2; di10 = 16.64; dP = 1.73; r0 = 6; m0 = 5 * 10^-2; n0 = 0.14; mG0 = 0.047; a0 = 0.14; kBarABeta0 = 7 * 10^-3; kBarNd = 10^-3; ki10 = 2.5 * 10^-6; kTBeta = 2.5 * 10^-7; kM = 0.047; kMHat = 0.047; kM1 = 0.03; kM2 = 0.017; kM1Hat = 0.04; kM2Hat = 0.007; kfI = 3.36 * 10^-10; kFO = 2.58 * 10^-11; kAO = 10^-7; kP = 6 * 10^-9; ktAlpha = 4 * 10^-5; (* System of Nonlinear PDEs *) pde ={D[aBetaI[x,y,t],t] == (lambdaBetaI* (1 +If[0<=t<=100,r0 * t / 100,r0]) - dABetaI * aBetaI[x,y,t]) * (n[x,y,t] / n0), D[aBeta0[x,y,t],t] == aBetaI[x,y,t] * Abs[D[n[x,y,t],t]] + lambdaN * (n[x,y,t] / n0) + lambdaA * (a[x,y,t] / a0) - (dABeta0MHat * (m1Hat[x,y,t] + myTheta * m2Hat[x,y,t]) + dABeta0M * (m1[x,y,t] + myTheta * m2[x,y,t])) * (aBeta0[x,y,t] / (aBeta0[x,y,t] + kBarABeta0)), D[myTau[x,y,t],t] == (lambdaTau0 + lambdaTau *If[0<=t<=100,r0 * t / 100,r0] - dTau * myTau[x,y,t]) * (n[x,y,t] / n0), D[fI[x,y,t],t] == (lambdaF * myTau[x,y,t] - dFI * fI[x,y,t])*(n[x,y,t] / n0), D[fO[x,y,t],t] == fI[x,y,t] * Abs[D[n[x,y,t],t]] - dFO * fO[x,y,t], D[n[x,y,t],t] == -dNF * (fI[x,y,t] / (fI[x,y,t] + kfI)) * n[x,y,t] - dNT * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) * (1 / (1 + myGamma * (i10[x,y,t] / ki10)))*n[x,y,t], D[a[x,y,t],t] == lambdaAABeta0 * aBeta0[x,y,t] + lambdaATAlpha * tAlpha[x,y,t] - dA * a[x,y,t], D[nD[x,y,t],t] == dNF * (fI[x,y,t] / (fI[x,y,t] + kfI)) * n[x,y,t] + dNT * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) * (1 / (1 + myGamma * (i10[x,y,t] / ki10)))*n[x,y,t] - dNdM * (m1[x,y,t] + m2[x,y,t]) * (nD[x,y,t] / (nD[x,y,t] + kBarNd)) - dNdMhat * (m1Hat[x,y,t] + m2Hat[x,y,t])*(nD[x,y,t] / (nD[x,y,t] + kBarNd)), D[aO[x,y,t],t] == lambdaAO * aBeta0[x,y,t] - dAO * aO[x,y,t] + DAO * Laplacian[aO[x,y,t],{x,y,t}], D[h[x,y,t],t] == lambdaH *nD[x,y,t] - dH * h[x,y,t] + DH * Laplacian[h[x,y,t],{x,y,t}], D[m1[x,y,t],t] == mG0 * (lambdaMF * (fO[x,y,t] / (fO[x,y,t] + kFO)) + lambdaMA * (aO[x,y,t] / (aO[x,y,t] + kAO))) * ((myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha))) /(myBeta *(tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) +(i10[x,y,t] / (i10[x,y,t] + ki10)))) - lambdaM1TBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1[x,y,t] - dM1 * m1[x,y,t] - Div[m1[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}], D[m2[x,y,t],t] == mG0 * (lambdaMF * (fO[x,y,t] / (fO[x,y,t] + kFO)) + lambdaMA * (aO[x,y,t] / (aO[x,y,t] + kAO))) * ((i10[x,y,t] / (i10[x,y,t] + kI10)) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + kI10)))) + lambdaM1TBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1[x,y,t] - dM2 * m2[x,y,t] - Div[m2[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}], D[m1Hat[x,y,t],t] == (alpha1 * (p[x,y,t] / (p[x,y,t] + kP))) * (m0 - ( m1Hat[x,y,t] + m2Hat[x,y,t])) * ((myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha))) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + ki10)))) - lambdaM1HatTBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1Hat[x,y,t] - dM1Hat * m1Hat[x,y,t] - Div[m1Hat[x,y,t] * Grad[aO[x,y,t],{x,y,t}],{x,y,t}], D[m2Hat[x,y,t],t] == (alpha1 * (p[x,y,t] / (p[x,y,t] + kP))) * (m0 - ( m1Hat[x,y,t] + m2Hat[x,y,t])) * ((i10[x,y,t] / (i10[x,y,t] + ki10)) /(myBeta * (tAlpha[x,y,t] / (tAlpha[x,y,t] + ktAlpha)) + (i10[x,y,t] / (i10[x,y,t] + ki10)))) + lambdaM1HatTBeta * (tBeta[x,y,t] / (tBeta[x,y,t] + kTBeta)) * m1Hat[x,y,t] - dM2Hat * m2Hat[x,y,t] - Div[m2Hat[x,y,t] * Grad[aO[x,y,t],{x,y,t}],{x,y,t}], D[tBeta[x,y,t],t] == lambdaTBetaM2 *m2[x,y,t] + lambdaTBetaM2Hat * m2Hat[x,y,t] - dTBeta * tBeta[x,y,t] + DTBeta * Laplacian[tBeta[x,y,t],{x,y,t}], D[tAlpha[x,y,t],t] == lambdaTAlphaM1 *m1[x,y,t] + lambdaTAlphaM1Hat * m1Hat[x,y,t] - dTAlpha * tAlpha[x,y,t] + DTAlpha * Laplacian[tAlpha[x,y,t],{x,y,t}], D[i10[x,y,t],t] == lambdai10M2 *m2[x,y,t] + lambdai10M2Hat * m2Hat[x,y,t] - di10 * i10[x,y,t] + Di10 * Laplacian[i10[x,y,t],{x,y,t}], D[p[x,y,t],t]== lambdaPA *a[x,y,t] + lambdaPM2 * m2[x,y,t] - dP * p[x,y,t] + DP * Laplacian[p[x,y,t],{x,y,t}]} ; (* Periodic Boundary Conditions *) bc = {aO[0,y,t] == aO[1,y,t], aO[x,0,t] == aO[x,1,t],h[0,y,t] ==h[1,y,t], h[x,0,t] == h[x,1,t], tBeta[0,y,t] == tBeta[1,y,t], tBeta[x,0,t] == tBeta[x,1,t], i10[0,y,t] == i10[1,y,t], i10[x,0,t] == i10[x,1,t], tAlpha[0,y,t] == tAlpha[1,y,t], tAlpha[x,0,t] == tAlpha[x,1,t], p[0,y,t] == p[1,y,t], p[x,0,t] == p[x,1,t]}; (* Initial Conditions *) ic = {aBetaI[x,y,0] == 10^(-6),aBeta0[x,y,0] == 10^(-8),myTau[x,y,0] == 1.37*10^(-10),fI[x,y,0] == 3.36*10^(-10),fO[x,y,0] == 3.36*10^(-11),n[x,y,0] == 0.14,a[x,y,0] == 0.14,m1[x,y,0] == 0.02,m2[x,y,0] == 0.02,m1Hat[x,y,0] == 0,m2Hat[x,y,0] == 0, nD[x,y,0] == 0,h[x,y,0] == 1.3*10^(-11),tBeta[x,y,0] == 10^(-6),tAlpha[x,y,0] == 2*10^{-5},i10[x,y,0] == 10^(-5),p[x,y,0] == 5*10^(-9),aO[x,y,0] == 0}; eqns = Flatten@{pde, bc, ic}; sols = NDSolve[eqns,{aBetaI,aBeta0,myTau,fI,fO,n,a,m1,m2,m1Hat,m2Hat,nD,h,tBeta,tAlpha,i10,p,aO},{x,0,1},{y,0,1},{t,0,3650}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 200}}];

I thought that I should be able to use the same method that is used here but I get the following error:

NDSolve: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method -> {"EquationSimplification"->"Residual"}.

Does anyone see any reason why I can't use the method that I'm trying to use? Also, should the variables that do not have laplacian terms only be a function of t? Currently, every variable is a function of x, y and t.

UPDATE: When I replace the last line with the following:

sols = NDSolve[eqns {aBetaI,aBeta0,myTau,fI,fO,n,a,m1,m2,m1Hat,m2Hat,nD,h,tBeta,tAlpha,i10,p,aO},{x,0,1},{y,0,1},{t,0,3650},Method ->{"EquationSimplification"->"Residual"}]

I get an error: enter image description here

Improve this question
$\endgroup$
9
  • $\begingroup$ Are you sure you wrote the system of equations 1-18 correctly? $\endgroup$ Commented Oct 16, 2019 at 16:13
  • $\begingroup$ I looked over them a couple of times to try to verify. I believe that I did but it's definitely a possibility that I have a typo. $\endgroup$ Commented Oct 16, 2019 at 17:05
  • $\begingroup$ Check first $\nabla.(M_1\nabla H)=(\nabla M_1)(\nabla H)+M_1\nabla^2 H$ (You have 4 typos there eq 11-14). And Abs(D[n[x,y,t],t]) should be Abs[D[n[x,y,t],t]] (You have 1 typo there eq 5)) $\endgroup$ Commented Oct 16, 2019 at 17:24
  • $\begingroup$ I fixed eq 5 by replacing the brackets and I fixed equations 11-14 by replacing the last term by Div[m1[x,y,t] * Grad[h[x,y,t],{x,y,t}],{x,y,t}] or the equivalent. I also used the Laplacian function for the 6 equations that have it instead of computing that by hand. Now, I'm getting the following error: NDSolve: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. $\endgroup$ Commented Oct 16, 2019 at 17:52
  • $\begingroup$ I also updated the code above. $\endgroup$ Commented Oct 16, 2019 at 17:55

1 Answer 1

Reset to default
2
$\begingroup$

After all the corrections, I also added artificial viscosity (but have not used it yet). The authors of the article did not indicate a solution method. They use average parameters. I used the parameters at point {.5,.5}. Therefore, there is no complete coincidence with Fig. 1 from the article.

(*Parameters*)DAO = 4.32*10^-2; aVis = 0 (*artificial viscosity *); DH = 8.11*10^-2; DTAlpha = 6.55*10^-2; DTBeta = 6.55*10^-2; Di10 = 6.04*10^-2; DP = 1.2*10^-1; lambdaBetaI = 9.51*10^-6; lambdaN = 8*10^-9; lambdaA = 8*10^-10; lambdaTau0 = 8.1*10^-11; lambdaTau = 1.35*10^-11; lambdaF = 1.662*10^-3; lambdaATAlpha = 1.54; lambdaAABeta0 = 1.793; lambdaAO = 5*10^-2; lambdaH = 3*10^-5; lambdaMF = 2*10^-2; lambdaMA = 2.3*10^-3; lambdaM1TBeta = 6*10^-3; lambdaM1HatTBeta = 6*10^-4; lambdaTBetaM2 = 1.5*10^-2; lambdaTBetaM2Hat = 1.5*10^-2; lambdaTAlphaM1 = 3*10^-2; lambdaTAlphaM1Hat = 3*10^-2; lambdai10M2 = 6.67*10^-3; lambdai10M2Hat = 6.67*10^-3; lambdaPA = 6.6*10^-8; lambdaPM2 = 1.32*10^-7; myTheta = 0.9; alpha1 = 5; myBeta = 10; myGamma = 1; dABetaI = 9.51; dABeta0M = 9.51; dABeta0M = 2*10^-3; dABeta0MHat = 10^-2; dTau = 0.277; dFI = 2.77*10^-3; dFO = 2.77*10^-4; dN = 1.9*10^-4; dNF = 3.4*10^-4; dNT = 1.7*10^-4; dNdM = 0.06; dNdMhat = 0.02; dA = 1.2*10^-3; dM1 = 0.015; dM2 = 0.015; dM1Hat = 0.015; dM2Hat = 0.015; dAO = 0.951; dH = 58.71; dTAlpha = 55.45; dTBeta = 3.33*10^2; di10 = 16.64; dP = 1.73; r0 = 6; m0 = 5*10^-2; n0 = 0.14; mG0 = 0.047; a0 = 0.14; kBarABeta0 = 7*10^-3; kBarNd = 10^-3; ki10 = 2.5*10^-6; kTBeta = 2.5*10^-7; kM = 0.047; kMHat = 0.047; kM1 = 0.03; kM2 = 0.017; kM1Hat = 0.04; kM2Hat = 0.007; kfI = 3.36*10^-10; kFO = 2.58*10^-11; kAO = 10^-7; kP = 6*10^-9; ktAlpha = 4*10^-5; fff[t_] := Piecewise[{{r*t/100, t <= 100}, {r, True}}] /. r -> r0; (*System of Nonlinear PDEs*) pde = {D[aBetaI[x, y, t], t] - ((lambdaBetaI*(1 + fff[t]) - dABetaI*aBetaI[x, y, t])*(n[x, y, t]/n0)) - aVis Laplacian[aBetaI[x, y, t], {x, y}], D[aBeta0[x, y, t], t] - (aBetaI[x, y, t]* Abs[-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] - dNT*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 + myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]] + lambdaN*(n[x, y, t]/n0) + lambdaA*(a[x, y, t]/ a0) - (dABeta0MHat*(m1Hat[x, y, t] + myTheta*m2Hat[x, y, t]) + dABeta0M*(m1[x, y, t] + myTheta*m2[x, y, t]))*(aBeta0[x, y, t]/(aBeta0[x, y, t] + kBarABeta0))) - aVis Laplacian[aBeta0[x, y, t], {x, y}], D[myTau[x, y, t], t] - ((lambdaTau0 + lambdaTau*fff[t] - dTau*myTau[x, y, t])*(n[x, y, t]/n0)), D[fI[x, y, t], t] - (lambdaF*myTau[x, y, t] - dFI*fI[x, y, t])*(n[x, y, t]/n0) - aVis Laplacian[fI[x, y, t], {x, y}], D[fO[x, y, t], t] - (fI[x, y, t]* Abs[-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] - dNT*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 + myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]] - dFO*fO[x, y, t]) - aVis Laplacian[fO[x, y, t], {x, y}], D[n[x, y, t], t] - (-dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] - dNT*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 + myGamma*(i10[x, y, t]/ki10)))*n[x, y, t]), D[a[x, y, t], t] - (lambdaAABeta0*aBeta0[x, y, t] + lambdaATAlpha*tAlpha[x, y, t] - dA*a[x, y, t]) - aVis Laplacian[a[x, y, t], {x, y}], D[nD[x, y, t], t] - (dNF*(fI[x, y, t]/(fI[x, y, t] + kfI))*n[x, y, t] + dNT*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha))*(1/(1 + myGamma*(i10[x, y, t]/ki10)))*n[x, y, t] - dNdM*(m1[x, y, t] + m2[x, y, t])*(nD[x, y, t]/(nD[x, y, t] + kBarNd)) - dNdMhat*(m1Hat[x, y, t] + m2Hat[x, y, t])*(nD[x, y, t]/(nD[x, y, t] + kBarNd))) - aVis Laplacian[nD[x, y, t], {x, y}], D[aO[x, y, t], t] - ((lambdaAO*aBeta0[x, y, t] - dAO*aO[x, y, t] + DAO*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))), D[h[x, y, t], t] - (lambdaH*nD[x, y, t] - dH*h[x, y, t] + DH*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}])), D[m1[x, y, t], t] - (mG0*(lambdaMF*(fO[x, y, t]/(fO[x, y, t] + kFO)) + lambdaMA*(aO[x, y, t]/(aO[x, y, t] + kAO)))*((myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)))/(myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y, t]/(i10[x, y, t] + ki10)))) - lambdaM1TBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))* m1[x, y, t] - dM1*(D[m1[x, y, t], x] D[h[x, y, t], x] + D[m1[x, y, t], y] D[h[x, y, t], y] + m1[x, y, t]*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}]))), D[m2[x, y, t], t] - (mG0*(lambdaMF*(fO[x, y, t]/(fO[x, y, t] + kFO)) + lambdaMA*(aO[x, y, t]/(aO[x, y, t] + kAO)))*((i10[x, y, t]/(i10[x, y, t] + ki10))/(myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y, t]/(i10[x, y, t] + ki10)))) + lambdaM1TBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))* m1[x, y, t] - dM2*(D[m2[x, y, t], x] D[h[x, y, t], x] + D[m2[x, y, t], y] D[h[x, y, t], y] + m2[x, y, t]*(D[h[x, y, t], {x, 2}] + D[h[x, y, t], {y, 2}]))), D[m1Hat[x, y, t], t] - ((alpha1*(p[x, y, t]/(p[x, y, t] + kP)))*(m0 - (m1Hat[x, y, t] + m2Hat[x, y, t]))*((myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)))/(myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y, t]/(i10[x, y, t] + ki10)))) - lambdaM1HatTBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))* m1Hat[x, y, t] - dM1Hat*(D[m1Hat[x, y, t], x] D[aO[x, y, t], x] + D[m1Hat[x, y, t], y] D[aO[x, y, t], y] + m1Hat[x, y, t]*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))), D[m2Hat[x, y, t], t] - ((alpha1*(p[x, y, t]/(p[x, y, t] + kP)))*(m0 - (m1Hat[x, y, t] + m2Hat[x, y, t]))*((i10[x, y, t]/(i10[x, y, t] + ki10))/(myBeta*(tAlpha[x, y, t]/(tAlpha[x, y, t] + ktAlpha)) + (i10[x, y, t]/(i10[x, y, t] + ki10)))) + lambdaM1HatTBeta*(tBeta[x, y, t]/(tBeta[x, y, t] + kTBeta))* m1Hat[x, y, t] - dM2Hat*(D[m2Hat[x, y, t], x] D[aO[x, y, t], x] + D[m2Hat[x, y, t], y] D[aO[x, y, t], y] + m2Hat[x, y, t]*(D[aO[x, y, t], {x, 2}] + D[aO[x, y, t], {y, 2}]))), D[tBeta[x, y, t], t] - (lambdaTBetaM2*m2[x, y, t] + lambdaTBetaM2Hat*m2Hat[x, y, t] - dTBeta*tBeta[x, y, t] + DTBeta*(D[tBeta[x, y, t], {x, 2}] + D[tBeta[x, y, t], {y, 2}])), D[tAlpha[x, y, t], t] - (lambdaTAlphaM1*m1[x, y, t] + lambdaTAlphaM1Hat*m1Hat[x, y, t] - dTAlpha*tAlpha[x, y, t] + DTAlpha*(D[tAlpha[x, y, t], {x, 2}] + D[tAlpha[x, y, t], {y, 2}])), D[i10[x, y, t], t] - (lambdai10M2*m2[x, y, t] + lambdai10M2Hat*m2Hat[x, y, t] - di10*i10[x, y, t] + Di10*(D[i10[x, y, t], {x, 2}] + D[i10[x, y, t], {y, 2}])), D[p[x, y, t], t] - (lambdaPA*a[x, y, t] + lambdaPM2*m2[x, y, t] - dP*p[x, y, t] + DP*(D[p[x, y, t], {x, 2}] + D[p[x, y, t], {y, 2}]))}; (*Periodic Boundary Conditions*) bc = {aO[0, y, t] == aO[1, y, t], aO[x, 0, t] == aO[x, 1, t], h[0, y, t] == h[1, y, t], h[x, 0, t] == h[x, 1, t], tBeta[0, y, t] == tBeta[1, y, t], tBeta[x, 0, t] == tBeta[x, 1, t], i10[0, y, t] == i10[1, y, t], i10[x, 0, t] == i10[x, 1, t], tAlpha[0, y, t] == tAlpha[1, y, t], tAlpha[x, 0, t] == tAlpha[x, 1, t], p[0, y, t] == p[1, y, t], p[x, 0, t] == p[x, 1, t]}; (*Initial Conditions*) ic = {aBetaI[x, y, 0] == 10^(-6), aBeta0[x, y, 0] == 10^(-8), myTau[x, y, 0] == 1.37*10^(-10), fI[x, y, 0] == 3.36*10^(-10), fO[x, y, 0] == 3.36*10^(-11), n[x, y, 0] == 0.14, a[x, y, 0] == 0.14, m1[x, y, 0] == 0.02, m2[x, y, 0] == 0.02, m1Hat[x, y, 0] == 0, m2Hat[x, y, 0] == 0, nD[x, y, 0] == 0, h[x, y, 0] == 1.3*10^(-11), tBeta[x, y, 0] == 10^(-6), tAlpha[x, y, 0] == 2*10^(-5), i10[x, y, 0] == 10^(-5), p[x, y, 0] == 5*10^(-9), aO[x, y, 0] == 0}; eqns = Flatten@{Table[pde[[i]] == 0, {i, 1, Length[pde]}], bc, ic}; var = {aBetaI, aBeta0, myTau, fI, fO, n, a, m1, m2, m1Hat, m2Hat, nD, h, tBeta, tAlpha, i10, p, aO}; Monitor[sols = NDSolve[eqns, var, {x, 0, 1}, {y, 0, 1}, {t, 0, 3650}, Method -> {"IndexReduction" -> Automatic, "EquationSimplification" -> "Residual", "PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 41, "MaxPoints" -> 41, "DifferenceOrder" -> 2}}}, EvaluationMonitor :> (monitor = Row[{"t=", CForm[t], " csol=", CForm[n[.5, .5, t]]}])], monitor];

Solution of the system of equations at point {x,y}={.5,.5}

point[i_] := With[{max = Max[Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]], min = Min[ Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]]}, Range[min, max, (-min + max)/4]] Table[Plot[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650}, PlotLabel -> var[[i]], PlotRange -> All, Ticks -> {{0, 2000, 3650}, point[i]}], {i, 1, Length[var]}]

Figure 1

Solution of the system of equations at point t=3650

Table[With[{max = Max[Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]], min = Min[ Table[var[[i]][.5, .5, t] /. sols[[1, i]], {t, 0, 3650, 50}]]}, Plot3D[var[[i]][x, y, 3650] /. sols[[1, i]], {x, 0, 1}, {y, 0, 1}, PlotLabel -> var[[i]], PlotRange -> {min, max}, Mesh -> None, ColorFunction -> Hue]], {i, 1, Length[var]}]

Figure2

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you so much for all of your help. I need the plots to match, so I will try to use the code you've provided to get matching plots. $\endgroup$ Commented Oct 18, 2019 at 13:26
  • $\begingroup$ @R'maniHaulcy You're welcome! $\endgroup$ Commented Oct 18, 2019 at 14:34

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.