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"}]
Abs(D[n[x,y,t],t])
should beAbs[D[n[x,y,t],t]]
(You have 1 typo there eq 5)) $\endgroup$ – Alex Trounev Oct 16 '19 at 17:24