Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

  1. 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 order n)
  2. 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.

share|improve this question

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.