0
$\begingroup$

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}]
$\endgroup$
  • $\begingroup$ Yes, it's wrong. Try executing e.g. y[x] /. x -> R == 0 and observe the output, and think about why you obtain this. $\endgroup$ – xzczd Feb 18 '19 at 14:28
  • $\begingroup$ 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 after mGB, but I cannot figure it out. Why the name of the functions in this list involve NDSolve? $\endgroup$ – MOON Feb 18 '19 at 14:55
  • $\begingroup$ Coding equations in this way makes them hard to check. I suggest using With like in e.g. this post to simplify the code. $\endgroup$ – xzczd Feb 18 '19 at 15:21
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Feb 18 '19 at 16:41
  • $\begingroup$ @MichaelE2 What is the difference between NDSolve and NDSolveValue? $\endgroup$ – MOON Feb 19 '19 at 14:48
1
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. Could you please let me know what the problem was? $\endgroup$ – MOON Feb 18 '19 at 21:30
  • $\begingroup$ I have removed all /.r->R options so they are not required in this task. $\endgroup$ – Alex Trounev Feb 18 '19 at 23:46

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.