I am solving a non linear partial differential equation with what I call free boundary conditions (solid mechanicists would know this as simply supported).
I realized that this boundary condition does induce a certain stiffness in my equation which is not so strong for, say, periodic boundary condition. How do I know this? The simulation runs longer for free boundaries with a stiff method such as LSODA
. This tells me that the time step is being made smaller and smaller to accommodate this rapidly changing time variation.
I also have various structures forming as obvious from the plot of my result. So I use the TensorProductGrid
spatial discretization method as I gather it would be useful to solve this equation with a preferential spatial grid treatment (denser grid in regions of necessity).
As I increase the DifferenceOrder
for my method/spatial discretization, I notice that I run into convergence issues NDSolve::ndcf
at a later stage. Viz., for a smaller DifferenceOrder
I have a convergence issue at a later time. This initially struck me as being a little strange as wouldn't a superior difference order help alleviate convergence? But then obviously, a greater difference order would catch the issue with convergence earlier!
Here is my code:
m=0.05;
\[Epsilon]=10^-6;
Bo=1/300;
\[Delta]=10^-7;
Bi=1;
K1=1;
wn=2;
(*R=Ra/(Pr S1);*)
R=0*1/10;
{xMin,xMax}=wn{0,79.5788};
k=0.0677*92.389/(xMax/wn);
TMax=10000*100;
n=1;
uSolpbc[t_,x_]=u[t,x]/.NDSolve[{\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x]\)\)+\!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((
\*SuperscriptBox[\(u[t, x]\), \(3\)]\
\*SubscriptBox[\(\[PartialD]\), \(x, x, x\)]u[t, x])\)\)-Bo \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((
\*SuperscriptBox[\(u[t, x]\), \(3\)]\
\*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\)+m \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((
\*SuperscriptBox[\((
\*FractionBox[\(u[t, x]\), \(K1 + Bi\ u[t, x]\)])\), \(2\)]\
\*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\) +\[Epsilon]/(Bi u[t,x] + K1) +\[Delta] u[t,x]^3/(Bi u[t,x] + K1)^3 D[D[u[t,x],x],x] + R \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((
\*SuperscriptBox[\(u[t, x]\), \(3\)]\
\*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\)== 0,
u[0,x]==1-0.05 (Cos[k*x]),
(*PERIODIC BOUNDARY CONDITIONS-START*)
(*u[t,xMin]== u[t,xMax],
Derivative[0,1]u[t,xMin]==Derivative[0,1]u[t,xMax],
Derivative[0,2]u[t,xMin]==Derivative[0,2]u[t,xMax],
Derivative[0,3]u[t,xMin]==Derivative[0,3]u[t,xMax]*)
(*PERIODIC BOUNDARY CONDITIONS-END. Uncomment to run periodic case*)
(*FREE BOUNDARY CONDITIONS-START*)
Derivative[0,1][u][t,xMin]==0,Derivative[0,1][u][t,xMax]==0,Derivative[0,3][u][t,xMin]==0,Derivative[0,3][u][t,xMax]==0
(*FREE BOUNDARY CONDITIONS-START*)
},
u,
{t,0,TMax},
{x,xMin,xMax},
(*PrecisionGoal->Automatic,*)
Method->{"MethodOfLines","Method"->"LSODA","TemporalVariable"->t,"SpatialDiscretization"->{"TensorProductGrid","MinPoints"->10000,"MaxPoints"->10000,"DifferenceOrder"->1}}][[1]]
And the plotting function:
evapfbc2 = Plot[
uSolpbc[6981*100, x], {x, xMin, xMax},
PlotStyle -> {Black, Thick, Dashed},
AxesLabel -> {"\[Lambda]", "h"},
BaseStyle -> {FontWeight -> "Plain", FontSize -> 18},
PlotRange -> {{0, xMax}, {0, 3.5}}
]
My questions are:
- Should I be increasing the difference order for better results? (I understand that with a difference order set to
n
, it allows the stiff solver "elbow room" up to and including ordern
) - Is this the wrong method to use and is there something I could do in mathematica that would make this better? I am assuming that I would need to use a stiff solver for this problem.
Note:
With a difference order of 1, this runs for about 10-15 seconds on an i7/8GB RAM system. With a difference order of 12, it runs for about 20 minutes and consumes about 3-4 gigs of RAM.