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]
```
Improve this question
$\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
    Commented Mar 14, 2020 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
    Commented Mar 14, 2020 at 14:40

2 Answers 2

Reset to default
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]
Improve this answer
$\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]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Sorry, @maeinss, i edited your code accidentally when trying to edit mine. Don't know how to undo. $\endgroup$
    – Akku14
    Commented Apr 15, 2020 at 6:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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