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.

Is the EventLocator option not compatible with LSODA on NDSolve. Below is what I tried to do to arrest my simulation if the time step went below a certain value (based on this suggestion).

The script file is also attached.

The error messages that I get are:

How should I go about it so that I can use EventLocator with LSODA. The reason I am doing this is because my equation turns stiff and the last half of the time steps are extremely small, slow and unnecessary for me.

NDSolve::moptx: Method option Method in {NDSolve`LSODA, 
MaxDifferenceOrder -> 12, Method -> {EventLocator, Event :>
If[stepsize < 10  , 0, 1]}}} is not one of {LinearSolveMethod,
MaxDifferenceOrder}.

NDSolve::initf: The initialization of the method NDSolve`LSODA failed.
#!/usr/local/bin/MathematicaScript -script

$HistoryLength=0;
    $pwf=$InputFileName;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0,EvapThickFilm,h,Bo,\[Epsilon],K1,\[Delta],Bi,m,r]
Eq0[h_,{Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_}]:=\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\)+Div[-h^3 Bo Grad[h]+h^3 Grad[Laplacian[h]]+(\[Delta] h^3)/(Bi h+K1)^3 Grad[h]+m (h/(K1+Bi h))^2 Grad[h]]+\[Epsilon]/(Bi h+K1) + (r)D[D[(h^2/(K1+Bi h)),x] h^3,x] ==0;
SetCoordinates[Cartesian[x,y,z]];
EvapThickFilm[Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_]:=Eq0[h[x,y,t],{Bo,\[Epsilon],K1,\[Delta],Bi,m,r}];
TraditionalForm[EvapThickFilm[Bo,\[Epsilon],K1,\[Delta],Bi,m,r]];


L=79.5788; TMax=12500*100;
Off[NDSolve::mxsst];
Clear[Kvar];
Kvar[t_]:=  Piecewise[{{1,t<=1},{2,t>1}}]
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol=h/.NDSolve[{
(*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)

EvapThickFilm[0,1*10^-6,1,0.001,1,2*0.025,0],
h[0,y,t]==h[L,y,t],
h[x,0,t]==h[x,L,t],
(*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)
h[x,y,0]==1+(-0.05 Cos[2\[Pi] x/L] -0.05 Sin[2 \[Pi] x/L])(Cos[2\[Pi] y/L])
},
h,
{x, 0, L},
{y,0, L},
{t, 0, TMax},


StepMonitor :> (
       laststep = thisstep; thisstep = t; 
       stepsize = thisstep - laststep;
       ),


Method->{"LSODA","MaxDifferenceOrder"->12, Method -> {"EventLocator", 
         "Event" :> (If[stepsize < 10^-4, 0, 1])} },
MaxStepFraction->1/50,
PrecisionGoal->3
][[1]]




{TMin,TRup}=InterpolatingFunctionDomain[hSol][[3]];
hSolGridData=InterpolatingFunctionValuesOnGrid[hSol];
hSolCoords=InterpolatingFunctionCoordinates[hSol];
finalStep=InterpolatingFunctionCoordinates[hSol][[3]];
(*rupture=NumberForm[N[TRup/100],6];*)
rupture=TRup/100;

$parameterfile=StringJoin[$pwf,".dat"];
Export[$parameterfile, {0, 100, 0, 0.0001, 35.1, 7.02, 0, 3, 1, 5, SetPrecision[rupture,5]}];
    $matfile=StringJoin[$pwf,".mat"];
    Export[$matfile,hSolGridData];
(*Exports time step data*)
$timefile=StringJoin[$pwf,"_time",".mat"];
Export[$timefile,InterpolatingFunctionCoordinates[hSol][[3]]];

hGrid = InterpolatingFunctionGrid[hSol];
{TMin,TRup}=InterpolatingFunctionDomain[hSol][[3]];
Length[hGrid];
{nX,nY,nT}=Drop[Dimensions[hGrid],-1];

fac=0.98;


$epsfile0=StringJoin[$pwf,"_0",".eps"];
$pngfile0=StringJoin[$pwf,"_0",".png"];
$epsfileRup=StringJoin[$pwf,"_TRup",".eps"];
$pngfileRup=StringJoin[$pwf,"_TRup",".png"];
ic=Plot3D[hSol[x,y,0*TRup],{x,0,L},{y, 0, L},

PlotRange->{{0,L},{0,L},{0,3.5}},
BaseStyle->{FontWeight->"Plain",FontSize->18},
PlotPoints->65,
ColorFunction->GrayLevel
];
Export[$epsfile0,ic,ImageSize->{350,350}];
    Export[$pngfile0,ic,ImageResolution->350];
rupProfile=Plot3D[hSol[x,y,fac*TRup],{x,0,L},{y, 0, L},

PlotRange->{{0,L},{0,L},{0,3.5}},
BaseStyle->{FontWeight->"Plain",FontSize->18},
PlotPoints->65,
ColorFunction->GrayLevel
];
Export[$epsfileRup,rupProfile,ImageSize->{350,350}];
    Export[$pngfileRup,rupProfile,ImageResolution->350];

dataxy=Import[$matfile];
    datat=Import[$timefile];


(*$epsfile0=StringJoin[$pwf,"_0_dft",".eps"];
$pngfile0=StringJoin[$pwf,"_0_dft",".png"];
$epsfileRup=StringJoin[$pwf,"_TRup_dft",".eps"];
$pngfileRup=StringJoin[$pwf,"_TRup_dft",".png"];
FDataFirst=Abs[Fourier[dataxy[[All,All,1]]]];
dftfirst=MatrixPlot[FDataFirst]*)
share|improve this question
3  
"EventLocator" is a controller method, while "LSODA" is an integration method. Witness for instance NDSolve[{y''[x] == x y[x], y[0] == 1, y'[0] == 0}, y, {x, -10, 0}, Method -> {"EventLocator", Method -> "LSODA", "Event" -> y'[x]}]. –  J. M. Sep 23 '12 at 13:46
    
@J.M. Hmmm.. I didn't think that through. –  drN Sep 23 '12 at 13:48
    
@J.M. I have a further question: when you don't include the difference order, you aren't able to limit your differentiation to that order? I realize that since my difference order itself was set to 12, it wouldn't matter if I ommitted it. –  drN Sep 23 '12 at 13:53
    
You can still use the "MaxDifferenceOrder" option; replace Method -> "LSODA" in the snippet above with Method -> {"LSODA", "MaxDifferenceOrder" -> 12}. Yes, options can be nested in this case. The numerical solution of differential equations is a complicated affair, after all. ;) –  J. M. Sep 23 '12 at 14:05
1  
BTW: NDSolve`LSODA // Options shows the defaults taken by the "LSODA" method value; as you can see, $12$ is indeed the maximum order taken by default. –  J. M. Sep 23 '12 at 14:06
show 4 more comments

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.