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
add comment

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.