0

I have this system of PDE that I want to solve. I reduced my original system of PDEs to the system below to capture the error (the actual code is after the image at the buttom):

enter image description here

Function f1 does not depend on r explicitly and it has only value at r = R. In real problem r is the radius of a sphere and f1 is some concentration at the surface of the sphere. If I delete the highlighted part there is no error. However, if I include the highlighted part I get this error and I do not know how to resolve it:

Transpose::nmtx: The first two levels of {f1,NDSolve`xs$1041} cannot be transposed.

The code is here:

R = 2;

(*equations*)
eqns = {(\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(f1[r, t]\)\) /. r -> R) == 
    f2[r, t] - 0.01*(f1[r, t] /. r -> R),
   \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(f2[r, t]\)\) == -1*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(f2[r, t]\)\))
   };

(*initial condition*)
intis = {f1[r, 0] == 0, 
   f2[r, 0] == 0.1};

(*boundary condition*)
bc = {
   f2[R, t] ==  0.1,
   (Evaluate[\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(f2[r, t]\)\)] /. 
      r -> R) == (f1[r, t] /. r -> R)
   };

 (*solving the equations*)
{f1, f2} = 
 usol = NDSolveValue[
   Flatten[{eqns, intis, bc}], {f1, f2}, {r, 0.0, R}, {t, 0, 0.01}]
MOON
  • 3,864
  • 23
  • 49
  • Yes, it's wrong. Try executing e.g. y[x] /. x -> R == 0 and observe the output, and think about why you obtain this. – xzczd Feb 18 '19 at 14:28
  • I corrected the assignment by using parenthesis (y[x] /.x->R) == 0 . However, now I get the error: Transpose::nmtx: The first two levels of {mT,mD,mB,mBG,NDSolvexs$1044,NDSolvexs$1043,NDSolvexs$1045} cannot be transposed.Here it seems something is wrong with the functions aftermGB`, but I cannot figure it out. Why the name of the functions in this list involve NDSolve? – MOON Feb 18 '19 at 14:55
  • Coding equations in this way makes them hard to check. I suggest using With like in e.g. this post to simplify the code. – xzczd Feb 18 '19 at 15:21
  • Try NDSolve instead of NDSolveValue. Maybe you can fix that first. -- Minor style remark: I think {f1, f2} = ... is a bad idea, because then what were variables in eqns etc. now have values, making those previously defined things invalid. It makes the code less robust, hard to play with and debug, imo. – Michael E2 Feb 18 '19 at 16:41
  • @MichaelE2 What is the difference between NDSolve and NDSolveValue? – MOON Feb 19 '19 at 14:48
  • In this case, the error messages printed. I suspect the problem indicated by NDSolve is the critical one and perhaps causes the others messages in NDSolveValue. I'm suggesting you look at the NDSolve message, which indicates the problem is ill-posed. I don't know how to fix it, so I'm suggesting you might come up with an appropriate change. Alex Trounev basically guesses at a possible way to change the problem, which does in fact create a new problem that NDSolve will solve. But I don't know if the new problem is the one you want to solve or not. – Michael E2 Feb 19 '19 at 19:51

1 Answers1

1

I debug the source code

R = 3.95; eps = 10^-3;
d2 = 0.03;
d3 = 11;
alpha1 = 0.2;
alpha2 = 0.12/60;
alpha3 = 1;
beta1 = 0.266;
beta2 = 0.28;
beta3 = 1;
gamma1 = 0.2667;
gamma2 = 0.35;
delta1 = 0.00297;
delta2 = 0.35;
 eq1 = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mT[r, theta, phi, 
      t]\)\) == (alpha1*mBG[r, theta, phi, t] + alpha2)*
             mD[r, theta, phi, t] - alpha3*mT[r, theta, phi, t] + 
           beta1*mBG[r, theta, phi, t]*cD[r, theta, phi, t] + 
           d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mT[r, theta, phi, 
           t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mT[r, theta, phi, t]\)\) + 
                 1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mT[r, theta, phi, 
           t]\)\)),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mD[r, theta, phi, 
      t]\)\) == -(alpha1*mBG[r, theta, phi, t] + alpha2)*
             mD[r, theta, phi, t] + alpha3*mT[r, theta, phi, t] + 
           beta2*cD[r, theta, phi, t] - beta3*mD[r, theta, phi, t] + 
           d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mD[r, theta, phi, 
           t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mD[r, theta, phi, t]\)\) + 
                 1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mD[r, theta, phi, 
           t]\)\)),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mB[r, theta, phi, t]\)\) == 
         gamma1*mT[r, theta, phi, t]*cB[r, theta, phi, t] - 
           gamma2*mB[r, theta, phi, t] - 
           delta1*mB[r, theta, phi, t]*cG[r, theta, phi, t] + 
           delta2*mBG[r, theta, phi, t] + 
     d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mB[r, theta, phi, 
           t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mB[r, theta, phi, t]\)\) + 
                 1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mB[r, theta, phi, 
           t]\)\)),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mBG[r, theta, phi, t]\)\) == 
         delta1*mB[r, theta, phi, t]*cG[r, theta, phi, t] - 
           delta2*mBG[r, theta, phi, t] + 
     d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mBG[r, theta, phi, 
           t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mBG[r, theta, phi, t]\)\) + 
                 1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mBG[r, theta, phi, 
           t]\)\))};
   eq2 = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cD[r, theta, phi, t]\)\) == 
        d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cD[r, theta, phi, t]\)\) + 
              2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cD[r, theta, phi, t]\)\) +
              1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cD[r, theta, phi, 
          t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cD[r, theta, phi, t]\)\) + 
              1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cD[r, theta, phi, 
          t]\)\)),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cB[r, theta, phi, t]\)\) == 
        d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cB[r, theta, phi, t]\)\) + 
              2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cB[r, theta, phi, t]\)\) +
              1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cB[r, theta, phi, 
          t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cB[r, theta, phi, t]\)\) + 
              1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cB[r, theta, phi, 
          t]\)\)),
      \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cG[r, theta, phi, t]\)\) == 
        d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cG[r, theta, phi, t]\)\) + 
              2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cG[r, theta, phi, t]\)\) +
              1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cG[r, theta, phi, 
          t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cG[r, theta, phi, t]\)\) + 
              1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cG[r, theta, phi, 
          t]\)\))
      };
intis = {mT[r, theta, phi, 0] == 0, mD[r, theta, phi, 0] == 0, 
   mB[r, theta, phi, 0] == 0, mBG[r, theta, phi, 0] == 0, 
   cD[r, theta, phi, 0] == 0.001, cB[r, theta, phi, 0] == 0.001, 
   cG[r, theta, phi, 0] == 0.001};
bc = {mT[r, theta, eps, t] == 0.001, 
      mT[r, theta, 2*Pi, t] == 0.001,
      mT[r, eps, phi, t] == 0.001,
      mT[r, Pi, phi, t] ==  0.001,
      mD[r, theta, eps, t] ==  0.001,
      mD[r, theta, 2*Pi, t] == 0.001,
      mD[r, eps, phi, t] == 0.001,
      mD[r, Pi, phi, t] ==  0.001,
      mB[r, theta, eps, t] ==  0.001, 
      mB[r, theta, 2*Pi, t] == 0.001,
      mB[r, eps, phi, t] == 0.001,
      mB[r, Pi, phi, t] ==  0.001,
      mBG[r, theta, eps, t] ==  0.001,
      mBG[r, theta, 2*Pi, t] == 0.001,
      mBG[r, eps, phi, t] == 0.001,
      mBG[r, Pi, phi, t] ==  0.001,
      cD[r, theta, eps, t] ==  0.001, 
      cD[r, theta, 2*Pi, t] == 0.001,
      cD[r, eps, phi, t] == 0.001,
      cD[r, Pi, phi, t] ==  0.001,
      cD[R, theta, phi, t] ==  0.001,
      d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cD[r, theta, phi, t]\)\) /. 
              r -> R) == (-(beta1*mBG[R, theta, phi, t] + beta2)*
              cD[R, theta, phi, t] + beta3*mD[R, theta, phi, t]),
      cB[r, theta, eps, t] ==  0.001,
      cB[r, theta, 2*Pi, t] == 0.001,
      cB[r, eps, phi, t] == 0.001,
      cB[r, Pi, phi, t] ==  0.001,
      cB[R, theta, phi, t] ==  0.001,
      d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cB[r, theta, phi, t]\)\) /. 
              r -> R) == (-gamma1*mT[R, theta, phi, t]*
       cB[R, theta, phi, t] +
             gamma2*mB[R, theta, phi, t]),
      cG[r, theta, eps, t] ==  0.001,
      cG[r, theta, 2*Pi, t] == 0.001,
      cG[r, eps, phi, t] == 0.001,
      cG[r, Pi, phi, t] ==  0.001,
      cG[R, theta, phi, t] ==  0.001,
      d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cG[r, theta, phi, t]\)\) /. 
              r -> R) == (-delta1*mB[R, theta, phi, t]*
       cG[R, theta, phi, t] +
             delta2*mBG[R, theta, phi, t])
      };

usol = NDSolveValue[
  Flatten[{eq1, eq2, intis, bc}], {mT, mD, mB, mBG, cD, cB, cG}, {r, 
   eps, R}, {theta, eps, Pi}, {phi, eps, 2*Pi}, {t, 0, 0.01}]

var = {mT, mD, mB, mBG, cD, cB, cG}; 
With[{r = R, t = .01}, 
  Table[Plot3D[
    usol[[i]][R, theta, phi, .01], {theta, eps, Pi}, {phi, eps, 2*Pi},
     PlotLabel -> var[[i]], PlotRange -> All, Mesh -> None, 
    ColorFunction -> Hue], {i, 1, Length[var]}]] // Quiet

fig1

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106