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$ – Alex Trounev Oct 16 '19 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$ – R'mani Haulcy Oct 16 '19 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$ – Alex Trounev Oct 16 '19 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$ – R'mani Haulcy Oct 16 '19 at 17:52
  • $\begingroup$ I also updated the code above. $\endgroup$ – R'mani Haulcy Oct 16 '19 at 17:55
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$ – R'mani Haulcy Oct 18 '19 at 13:26
  • $\begingroup$ @R'maniHaulcy You're welcome! $\endgroup$ – Alex Trounev Oct 18 '19 at 14:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.