0
$\begingroup$

I'm new to Mathematica and have a hard time solving a stiff PDE using NDSolve after searching around for different stiff system methods. Below is the PDE and so far I have been getting the errors of NDSolve::ndsz: At t == xxx, step size is effectively zero; singularity or stiff system suspected. It's been quite frustrating and any feedback is deeply appreciated.

(*parameter setup*)
pm = 0.64; 
p0 = 0.5;
tmax = (pm - p0)/pm;
Tanhm1[t_, t0_] := Tanh[500*(t - t0)] + 1;
T[t_] := 293 + 20*Tanhm1[t, 0.1*tmax];
R = 225*^-9; k = 1.381*^-23;
mu[t_] := (2.414*^-5)*10^(247.8/(T[t] - 140)); 
g = 9.8; rp = 1050; rl = 1000;
D0[t_] := k*T[t]/(6*Pi*mu[t]*R);
U0[t_] := 2*R^2*g*(rp - rl)/(9*mu[t]);
H0 = 28*^-6; pinf = 0;
A1 = 8.07131; B1 = 1730.63; C1 = 233.426;
pint[t_] := 10^(A1 - B1/(C1 + T[t]));
kg = 20*^-16;
delta[t_] := kg*H0/(rl *D0[t]);
E0[t_] := delta[t]*(pint[t] - pinf);
Pe[t_] := E0[t]*H0/D0[t];
Ns[t_] := U0[t]/E0[t];
K[p_] := (1 - p)^6.55;
Z[p_] := 1/(pm - p);
DP [p_] := D[p*Z[p], p];

(*PDE*)
pde = {D[p[t, x], t] + (x/(1 - t))*D[p[t, x], x] == 
    Ns[t]/(1 - t)*D[K[p[t, x]]*p[t, x], x] + 
     1/(Pe[t]*(1 - t)^2)*D[K[p[t, x]]*DP[p[t, x]]*D[p[t, x], x], x]};
bc = {Ns[t]*Pe[t]*p[t, 0]*(1 - t)*(1 - Exp[-1000*t]) + 
     DP[p[t, 0]]*Derivative[0, 1][p][t, 0] == 0, 
   Ns[t]*Pe[t]*K[p[t, 1]]*p[t, 1]*(1 - t)*(1 - Exp[-1000*t]) + 
     K[p[t, 1]]*DP[p[t, 1]]*Derivative[0, 1][p][t, 1] == 
    Pe[t]*(1 - t)*p[t, 1]*(1 - Exp[-1000*t]),
   p[0, x] == p0};
sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MinPoints" -> 41, "MaxPoints" -> 101, 
      "DifferenceOrder" -> "Pseudospectral"}}]
Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, PlotPoints -> 100]
```
$\endgroup$
2
  • 1
    $\begingroup$ Removing "DifferenceOrder" -> "Pseudospectral" seems to help. Also, add the option PlotRange -> All in Plot3D to see the entire solution. $\endgroup$ – bbgodfrey Mar 14 '20 at 14:38
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Mar 14 '20 at 14:40
1
$\begingroup$

Removing all the options to NDSolve, it works fine in version 8.0. Simpliy do

sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1}]

Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, 
  PlotPoints -> 100, PlotRange -> All]
$\endgroup$
0
$\begingroup$

Can you use the StiffnessSwithching method?

 sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1}, 
       Method -> {"StiffnessSwitching"}]

 Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, 
       PlotPoints -> 100, PlotRange->All]
$\endgroup$
1
  • $\begingroup$ Sorry, @maeinss, i edited your code accidentally when trying to edit mine. Don't know how to undo. $\endgroup$ – Akku14 Apr 15 '20 at 6:43

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.