(*********************************************************************** Mathematica-Compatible Notebook This notebook can be used on any computer system with Mathematica 4.0, MathReader 4.0, or any compatible application. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. ***********************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 28961, 851]*) (*NotebookOutlinePosition[ 29849, 883]*) (* CellTagsIndexPosition[ 29773, 877]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["SplitStep method package", "Title"], Cell["\<\ This notebook is an annotated version of the package SplitStep.m.\ \>", "Text"], Cell["\<\ (* Package for using the split-step method for solving the cubic nonlinear Schrodinger equation D[u[x,t],t] + D[u[x,t],x,x] + 2 Abs[u[x,t]]^2 u[x,t] == v(x) u[x,t] *)\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["BeginPackage[\"SplitStep`\"];", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "What this does is sets up the context path so that any new symbols created \ (until a matching End[]) are created in the SplitStep` context. When I load \ the package, SplitStep` will be put onto the context path in the current ", StyleBox["Mathematica", FontSlant->"Italic"], " session, so that you can access these symbols without any extra typing. \ Thus, any new symbol which you want to have effectively exported to the ", StyleBox["Mathematica", FontSlant->"Italic"], " session at large needs to be created in this block. A symbol is created \ at the time it is first parsed, so something as simple as a list with the \ symbols you wanted to export right at the beginning would be adequate. What \ most people do (and what you will typically find in ", StyleBox["Mathematica", FontSlant->"Italic"], "'s standard packages is to in effect declare those symbols by giving their \ usage messages as I have done below." }], "Text"], Cell["\<\ SplitStepSolve::usage = \"SplitStepSolve[initialcondition, {x, xmin, xmax, \ n}, {{tmin, tmax, toutput}, deltat}] computes the solution of the cubic \ nonlinear Schrodinger equation on the spatial interval xmin to xmax from tmin \ to tmax using time steps of deltat and producing output at intervals toutput.\ \"; Options[SplitStepSolve] = {Potential->0, Order->4,Output->AbsValue, \ WorkingPrecision->$MachinePrecision, AccuracyGoal->Automatic, \ PrecisionGoal->Automatic,StoppingCondition->None, MovingWindow->False, \ AbsorbingBoundary->False}; Discretize::usage = \"Discretize[f,{x,xmin,xmax}, n] will produce a list of n \ values of f[x] from f[xmin] up to, but not including f[xmax].\"; Mass::usage = \"Mass[Abs[u], L] gives the mass analogue for the solution \ represented by the vector u over the length L.\"; Momentum::usage = \"Momentum[u,L] gives the momentum analogue for the \ solution represented by the vector u over the length L.\"; Centroid::usage = \"Centroid[Abs[u], grid] gives the position of the centroid \ of the discretized solution u on the given grid.\"; Variance::usage = \"Variance[Abs[u],grid] gives the variance of the solution \ represented by u on the given grid.\"; Value::usage = \"Value is an output option for SplitStepSolve which gives the \ solution vector.\"; AbsValue::usage = \"AbsValue is an output option for SplitStepSolve which \ gives the modulus of the solution vector.\"; MaxValue::usage = \"MaxValue is an output option for SplitStepSolve which \ gives the maximum modulus of the solution vector.\"; MaxPosition::usage = \"MaxPosition is an output option for SplitStepSolve \ which gives the position of maximum modulus of the solution vector.\"; Grid::usage = \"MaxPosition is an output option for SplitStepSolve which \ gives the current grid. This is useful to have when MovingWindow->True.\"; Time::usage = \"Refers to the output time at which values are given.\"; SplitStepSolve::precg = NIntegrate::precg; SplitStepSolve::accg = NIntegrate::accg; SplitStepSolve::timest = \"The time step `1` should be a positive real \ number.\"; SplitStepSolve::fixts = \"The AccuracyGoal and PrecisionGoal options have no \ effect with fixed time steps.\"; SplitStepSolve::output = \"`1` is not a valid output option for \ SplitStepSolve.\"; SplitStepSolve::ctrmov = \"You have chosen a moving window. Without either \ the centroid or the grid in the output, you will not be able to know the real \ position of the solution.\"; SplitStepSolve::absorb = \"The value of the option AbsorbingBoundary->`1` \ should be False, True, or a list of 2 numbers specifying the size and extent \ of the absorbing boundary condition.\"; SplitStepSolve::order = \"The value of the option Order->`1` should be either \ 1 or and even integer giving the temporal convergence order.\"; SplitStepSolve::stopt = \"The computation was not completed because the value \ of the option StoppingCondition->`1` evaluated to neither True nor False at \ `2`.\";\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ Certainly I do not want to export for general use all of the symbols which I \ define in the package. The use of a subcontext, here Private` (this is not a \ keyword) keeps all of the symbols used in definitions of the symbols \ \"declared\" in the context SplitStep`Private` which will not be seen unless \ their full name is given.\ \>", "Text"], Cell[CellGroupData[{ Cell["Begin[\"`Private`\"]", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[BoxData[ \("SplitStep`Private`"\)], "Output"] }, Open ]], Cell[TextData[{ "This is the definition for the linear part of the operator. A key feature \ of this is the use of a \"cache\" for the vectors which represent the \ differential operator in spectral space, which ", StyleBox["Mathematica", FontSlant->"Italic"], " makes particularly easy to implement." }], "Text"], Cell["\<\ Linear[dt_, L_][u_List] := InverseFourier[Fourier[u] CachedD2Vector[dt, \ Length[u], L]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ ExpFourierD2Vector is the basic function which generates the vector which \ represents the action of the differential operator in spectral space over a \ time interval dt with a spatial interval of extent L.\ \>", "Text"], Cell["\<\ ExpFourierD2Vector[dt_, n_, L_] := Block[{t, s = -I dt (2 Pi/L)^2}, \tt = Table[Exp[s j^2], {j, 1, n/2 - 1}]; \tJoin[{1}, t, {Exp[s (n/2)^2]}, Reverse[t]] ]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "Since the generator may be used many times (especially with automatic time \ step control), it is important to have it fast when possible. The best way \ to ensure this is to compile it as a part of loading the package. A ", ButtonBox["speed comparison", ButtonData:>{"CompiledSpeed.nb", None}, ButtonStyle->"Hyperlink"], " shows that this is well worth the effort." }], "Text"], Cell["\<\ CompiledExpFourierD2Vector = Compile[{{dt, _Real}, {n, _Integer}, {L, _Real}}, \tBlock[{t, s = -I dt (2 N[Pi]/L)^2}, \t\tt = Table[Exp[s j^2], {j, 1, n/2 - 1}]; \t\tJoin[{1.}, t, {Exp[s (n/2)^2]}, Reverse[t]] \t] ];\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "While having the compiled version is a lot faster, in a typical \ calculation, these vectors are used repeatedly. While it would be possible \ to save lists of these vectors and pass the lists to the functions called, \ that would be cumbersome and require more sophisticated preprocessing prior \ to each change in the time step. A cache does take up more memory (with some \ extra work, this could be limited to a specific amount), but saves doing the \ calculations over and over again. The caching is done by having the general \ definition of the function include a specific assignment for the given \ arguments. When using the rules given to define a function, ", StyleBox["Mathematica", FontSlant->"Italic"], " tries the specific before the general, so in subsequent calls with the \ same arguments, the cached value will be used. " }], "Text"], Cell[BoxData[ \(\(ClearD2Cache\ := \ \((\[IndentingNewLine]Clear[ CachedD2Vector]; \[IndentingNewLine]\((CachedD2Vector[dt_, n_, L_] := \((CachedD2Vector[dt, n, L] = If[MachineNumberQ[dt] || MachineNumberQ[L], CompiledExpFourierD2Vector[dt, n, L], ExpFourierD2Vector[dt, n, L]])\))\)\[IndentingNewLine])\); \)\)], "Input"], Cell["\<\ The first derivative is used only for processing the solution, so no attempt \ has been made to optimize it.\ \>", "Text"], Cell["\<\ FourierDVector[ord_, n_, L_] := Block[{t, s}, \tIf[ord === 0, Return[Table[1,{n}]]]; \ts = (2 Pi I/L)^ord; \tt = Table[s j^ord, {j, 1, n/2 - 1}]; \tJoin[{0}, t (-1)^ord, {s (n/2)^ord}, Reverse[t]] ]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ Now the nonlinear part of the operator. This has been written in a general \ way so it can be applied to various types of nonlinearities and potentials. \ This is based on a simple linearization: i.e. as if v[u] is constant.\ \>", "Text"], Cell["\<\ (* \tThe potential is assumed to be a nonlinear function of v \tfor the basic cubic NLSE, use, e.g. \tv = Function[2 # Conjugate[#]] *) Nonlinear[dt_, v_][u_List] := u Exp[I dt v[u]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ The basic first order split-step method involves simply alternating \ application of the linear and nonlinear parts of the approximate operator.\ \>", "Text"], Cell["\<\ SplitStep[1, dt_, v_, L_][u_] := Linear[dt, L][Nonlinear[dt, v][u]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ If you divide the linear steps at either end in half, you get second order in \ time. Note that with repeated application, the two half linear steps could be \ combined into one full one because of linearity.\ \>", "Text"], Cell["\<\ SplitStep[2, dt_, v_, L_][u_] := \tLinear[dt/2, L][Nonlinear[dt, v][Linear[dt/2,L][u]]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "Now we can use the method of Yoshida to go to arbitrary order. This will, \ in general work for this type of splitting for any equation which is \ time-reversible. Yoshida proposes alternate methods which save steps for \ order 6 and 8. I believe that with some effort, this could be generalized \ using ", StyleBox["Mathematica", FontSlant->"Italic"], " to generate step saving coefficients for arbitrary order, but I have not \ done this here because I have found in practice that the fourth order method \ typically gives the best efficiency for machine numbers." }], "Text"], Cell["\<\ (* \tThis is the method of Yoshida. \tDepends on the time-reversability of the NLSE *) SplitStep[n_?EvenQ, dt_, p___][u_] := Block[{z0, z1, alpha, beta, k = 1/(n - 1)}, \tz1 = 1/(2 - 2^k); \tz0 = - z1 2^k; \tSplitStep[n - 2, z1 dt, p][SplitStep[n - 2, z0 dt, p][SplitStep[n - 2, z1 \ dt, p][u]]]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "The method produced above will not be in an ideal form to use numerically. \ The most important thing that can easily be done it to combine two adjacent \ linear steps into one. The ", ButtonBox["simplified methods", ButtonData:>{"SplitStepSimplify.nb", None}, ButtonStyle->"Hyperlink"], " save several FFTs. It is worth noting that I have written the rules for \ SplitStep and SplitStepSimplify in a general way. The reason for doing this \ is that the splitting method can be applied to a variety of equations. An ", ButtonBox["example", ButtonData:>{"SplitStepSimplify.nb", "ODExample"}, ButtonStyle->"Hyperlink"], " shows that with a simple redefinition of Linear and Nonlinear, you can \ use the arbitrary order methods in a more general way." }], "Text"], Cell["\<\ (* \tTake the SplitStep method for order ord and simplify it for \tactual numerical use. *) SimplifySplitStep[ord_?EvenQ, {dt_, p___}, prec___] := Block[{Linear, Nonlinear, u}, \t(* This MUST be Block here to work. Since the methods are defined \t\tin terms of Linear (and Nonlinear), using Block allows you to \t\tgive a temporary definition for it without having that affect \t\tthe main one you want to be used in numerical evaluation (or \t\thaving the definition for numerical evaluation evaluate here.) \t*) \tModule[{new}, (* Use Module for local variables *) \t\t(* We make a (temporary) definition for Linear which will combine \t\t two Linear steps into a single one. *) \t\tLinear[dt1_, pl___][Linear[dt2_, pl___][x__]] := \t\t\tLinear[Simplify[dt1 + dt2], pl][x]; \t\tnew = SplitStep[ord, dt, p][Slot[1]]; \t\tnew = N[new, prec]; \t\t(* Now we want to make a Function that we can call to invoke a step \t\t of the simplified composite method. Using Function here is fast \t\t and it prevents evaluation until time of use. *) \t\targs = {dt, p}; \t\tFunction[Evaluate[args], Evaluate[Function[Evaluate[new]]]] \t] ]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ This is an exported function, which is just a convenience for using the \ package.\ \>", "Text"], Cell["\<\ Discretize[f_, {x_,xmin_, xmax_}, n_] := Block[{fx = Function[x,f]}, Table[fx[xmin + i(xmax - xmin)/n],{i,0,n-1}] ]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "Here is the definition of the main function of the package. I could spend \ hours going into details of all of the features I have packed into this \ function, so I will focus on just a couple of aspects, the automatic step \ size control and the construction of the output. A good example to \ understand the automatic step control is the ", ButtonBox["collision", ButtonData:>{"Collision.nb", None}, ButtonStyle->"Hyperlink"], " of two solitons. Investigating the action of ", ButtonBox["absorbing boundary condition potentials", ButtonData:>{"Absorbing.nb", None}, ButtonStyle->"Hyperlink"], " gives a good opportunity to view the construction of the output." }], "Text"], Cell["\<\ SplitStepSolve[InitialCondition_, {x_, xmin_, xmax_, n_}, {{tbegin_, tend_, \ T_:1},dt_}, opts___] := SplitStepSolve[Discretize[InitialCondition, \ {x,xmin,xmax}, n], {x, xmin, xmax}, {{tbegin, tend, T}, dt}, opts]; SplitStepSolve[InitialCondition_List, {x_, xmin_, xmax_}, {{tbegin_, tend_, \ T_:1},dtin_}, opts___] := Block[{n = Length[InitialCondition], \t\tL = xmax - xmin, \t\tpot, rv, v, ndt, order, out, wp, alpha, beta, stoptest, absorb, \t\tmoving,grid,deltax,position, last, \t\t$MaxPrecision, $MinPrecision}, \twp = WorkingPrecision /. {opts} /. Options[SplitStepSolve]; \tIf[Or[dtin === Automatic, And[ListQ[dtin], Length[dtin] > 0, dtin[[1]] === \ Automatic]], \t\tdt = Automatic; \t\tdts = N[If[And[ListQ[dtin], Length[dtin] > 1],dtin[[2]], T/1000], wp]; \t\tIf[Not[NumberQ[dts]], \t\t\tMessage[SplitStepSolve::timest, dts]; \t\t\tReturn[$Failed] \t\t], \t(* else *) \t\tdt = N[dtin, wp]; \t\tIf[Not[NumberQ[dts]], \t\t\tMessage[SplitStepSolve::timest, dts]; \t\t\tReturn[$Failed] \t\t] \t]; \t(* \t\tProcess the Accuracy and Precision goal for Automatic time step control \t*) \tIf[dt === Automatic, \t\tatol = AutoAGPG[wp,AccuracyGoal /. {opts} /. Options[SplitStepSolve], \ Accuracy]; \t\tIf[atol === $Failed, Return[$Failed]]; \t\trtol = AutoAGPG[wp, PrecisionGoal /. {opts} /. Options[SplitStepSolve], \ Precision]; \t\tIf[rtol === $Failed, Return[$Failed]], \t(* else *) \t\tIf[MemberQ[{opts}, AccuracyGoal | PrecisionGoal,Infinity], \t\t\tMessage[SplitStepSolve::fixts]] \t];\t\t\t \t(* \t\tThe Potential. Unless it is already a list, we discretize it \t*) \tpot = Potential /. {opts} /. Options[SplitStepSolve]; \tIf[!ListQ[pot], \t v = Discretize[pot,{x,xmin,xmax},n], \t v = pot]; \t(* \t\tThe output to keep track of \t*) \tout = Output /. {opts} /. Options[SplitStepSolve]; \tIf[!ListQ[out],out = {out}]; \tIf[Catch[Scan[ \t\t\tIf[!MatchQ[#,Value | AbsValue | MaxValue | MaxPosition | Mass | \ Momentum | Centroid | Variance | Grid | Plot | {Plot, \ ___}],Message[SplitStepSolve::output, #]; Throw[True]]&, \t\t\tout];False], \t\tReturn[$Failed]]; \tmoving = TrueQ[MovingWindow /. {opts} /. Options[SplitStepSolve]]; \tIf[moving, \t If[Not[MemberQ[out, Centroid | Grid]], \t\tMessage[SplitStepSolve::ctrmov]]; \t If[ListQ[pot], \t\tppot = True; \t\tMessage[SplitStepSolve::periodicv], \t\tpder = pot; \t\tppot = Catch[Do[ \t\t If[(pder /. x->xmax) != (pder /. x->xmin),Throw[False]]; \t\t pder = D[pder, x], \t\t {5}]; True] \t ]; \t]; \t(* \t\tThe is a simple default for AbsorbingBoundary->True which \t\tworks pretty well most of the time. \t*) \tabsorb = AbsorbingBoundary /. {opts} /. Options[SplitStepSolve]; \tIf[absorb =!= False, \t If[TrueQ[absorb], absorb = {1,1/4}]; \t If[Not[ListQ[absorb] && (Length[absorb] == 2)], \t\t\tMessage[SplitStepSolve::absorb, absorb]; \t\t\tReturn[$Failed]; \t ]; \t absorb = Discretize[AbsorbingBoundaryConditionPotential[absorb[[1]], \ {xmin,xmax}, {xmin, xmax} + L absorb[[2]] {1,-1}][x],{x,xmin,xmax},n]; \t]; \t(* \t\tThe order should be 1 or even \t*) \torder = Order /. {opts} /. Options[SplitStepSolve]; \tIf[Not[Or[order === 1, EvenQ[order]]], \t\tMessage[SplitStepSolve::order]; \t\tReturn[$Failed]; \t]; \tIf[wp === $MachinePrecision, np = Sequence[], np = wp]; \t(* \t\tDerive the integration method \t*) \tsplits = SimplifySplitStep[order, {hsym, vsym, Lsym},np]; \tgrid = Discretize[x,{x,xmin,xmax},n]; \toriggrid = grid; \tdeltax = (xmax - xmin)/n; \tposition = Round[(xmax + xmin)/(2*deltax)]; \tIf[wp === $MachinePrecision, \t\t(* \t\t\tWith machine precision, we set up complex PackedArrays \t\t\tthis ensures top speed. \t\t*) \t\tgrid = Developer`ToPackedArray[grid, Real]; \t\tulast = Developer`ToPackedArray[#, Complex]& /@ InitialCondition; \t\tv = Developer`ToPackedArray[v, Complex]; \t\tIf[absorb =!= False, \t\t\tabsorb = Developer`ToPackedArray[absorb, Complex]; \t\t\tav = v + absorb, \t\t(* else *) \t\t\tav = v], \t(* else \t\t\totherwise, we numericize to the desired precision \t\t*) \t\tgrid = N[grid, wp]; \t\tulast = N[InitialCondition, wp]; \t\tv = N[v, wp]; \t\tIf[absorb =!= False, \t\t\tabsorb = N[absorb, wp]; \t\t\tav = N[v + absorb + v, wp], \t\t(*else*) \t\t\tav = v]; \t\t(* \t\t\tFix the precision. These have been localized inside Block, so \t\t\tthis will not affect your session. For iterative processes \t\t\tit is often a good idea to fix the precision. \t\t*) \t\t$MaxPrecision = $MinPrecision = wp \t]; \tstoptest = StoppingCondition /. {opts} /. Options[SplitStepSolve]; \tIf[stoptest === None, stoptest = Function[True]]; \tglast = grid; \tvext = v; \ttt = tbegin; \tIf[dt === Automatic, \t\th = dts, \t(* else *) \t\tClearD2Cache; \t\tndt = Ceiling[T/dt]; \t\th = N[T/ndt, np]]; \tCatch[ \t\tdatabag = Internal`Bag[]; \t\tstepbag = Internal`Bag[]; \t\tInternal`StuffBag[databag,StepOutput[tbegin,ulast,out, grid]]; \t\tDo[ \t\t last = Internal`BagPart[databag, -1]; \t\t\tstopd = stoptest[last]; \t\t\tIf[stopd =!= False, \t\t\t\tIf[stopd =!= True, \t\t\t\t\tMessage[SplitStepSolve::stopt, stoptest,Rule[Time, Time /. last]]]; \t\t\t\tBreak[]]; \t\t\tIf[moving, \t\t\t\tcent = Centroid /. last; \t\t\t shift = Round[(cent/deltax) - position]; \t\t\t\tIf[shift != 0, \t\t\t\t\tposition += shift; \t\t\t\t\tgrid = origgrid + position*deltax; \t\t\t\t\tulast = RotateLeft[ulast, shift]; \t\t\t\t\tIf[ppot, \t\t\t\t\t\tv = RotateLeft[v, shift], \t\t\t\t\t(* else *) \t\t\t\t\t\tv = Discretize[pot,{x,xmin + position*deltax, \ xmin+position*deltax}] \t\t\t\t\t]; \t\t\t\t\tav = If[absorb === False,v, v+absorb]; \t\t\t\t]; \t\t\t]; \t\t\tIf[Developer`ZeroQ[Max[Abs[av]]], \t\t\t\tvfun = Function[2 # Conjugate[#]], \t\t\t\tvfun = Function[2 # Conjugate[#] - av] \t\t\t]; \t\t\tIf[dt === Automatic, \t\t\t\tClearD2Cache; \t\t\t\tWhile[Not[Or[tt == t, tt > t]], \t\t\t\t\tIf[tt + 2 h >= t, \t\t\t\t\tIf[Or[tt + h == t, tt + h > t], \t\t\t\t\t\th = t - tt, \t\t\t\t\t\t(* Avoid a very small step *) \t\t\t\t\t\th = (t - tt)/2 \t\t\t\t\t\t]; \t\t\t\t\t]; \t\t\t\t\t{h,hused,ulast} = DoubleStep[Function[{hh},Evaluate[splits[hh,vfun, \ L]]], h, ulast, order, {atol, rtol}]; \t\t\t\t\ttt += hused;\t\t\t\t \t\t\t\t\tInternal`StuffBag[stepbag,{tt,hused}]; \t\t\t\t], \t\t\t(* else: fixed step size *) \t\t\t\tulast = Nest[splits[h, vfun, L], ulast, ndt]; \t\t\t\ttt += h*ndt; \t\t\t\tInternal`StuffBag[stepbag,{tt,h}]; \t\t\t]; \t\t\tglast = grid; \t\t\tInternal`StuffBag[databag, StepOutput[t,ulast,out, grid]], \t\t\t{t, tbegin + T, tend, T}]; \t]; \tClearD2Cache; \tInternal`BagPart[databag, All] ]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell[TextData[{ "This is the module which handles automatic step control. As arguments, it \ eccepts a method, a proposed step size, the current solution, the method \ order, and the local tolerance criteria. The basic idea is to compute the \ solution using both one step of size h and two steps of size h/2 and to use \ Richardson extrapolation to make an estimate of the error. Suppose the method \ has error of order p. Then, if ", Cell[BoxData[ \(TraditionalForm\`A\_h\)]], "refers to the approximation with step size h, and S the actual solution, \ we have asymptotically\n", Cell[BoxData[ \(TraditionalForm\`\(\(A\_h = S\ + \ c\ h\^\(p + 1\)\)\(\ \)\)\)]], " and ", Cell[BoxData[ \(TraditionalForm\`A\_\(h/2\) = S\ + 2\ c\ \((h/2)\)\^\(p + 1\)\)]], "\n(p + 1 because this is the local error: combining steps gives one order \ less, namely p. Solving these equations for c, the error estimate is\n2 ", Cell[BoxData[ \(TraditionalForm \`\(c(h/2)\)\^p \[TildeFullEqual] \ \((A\_h - A\_\(h/2\))\)/\((2\^p - 1)\)\)]], ".\nIn a similar way, an \"optimal\" next step size h can also be computed. \ In practice we multiply this by a safety factor because it is relatively \ expensive to reject a step." }], "Text"], Cell["\<\ The scaled norm here is failrly typical for DE solvers. The balance of \ absolute (AccuracyGoal) and relative (PrecisionGoal) tolerance is the same as \ for the built in NDSolve command.\ \>", "Text"], Cell["\<\ DoubleStep[split_, htry_, u_, ord_, {atol_, rtol_}] := Block[{h = 2 htry, half, one, two, errest = Infinity, count = 0, fac = 1}, half = split[h/2][u]; While[errest > 1, \t\tcount++; \t\th /= 2; \t\tone = half; \t\thalf = split[h/2][u]; \t\ttwo = split[h/2][half]; \t\tw = atol + rtol*Map[Max, Transpose[{Abs[u], Abs[two]}]]; \t\tdiff = Max[Abs[(two - one)/w]]; \t\terrest = diff/(2^ord - 1); ]; (* Success: increase h if we didn't have to reduce *) If[count < 2, \tIf[Developer`ZeroQ[errest], \t fac = 2, \t fac = Min[2,(9/10)*(1/errest)^(1/(ord + 1))] \t] ]; {h*fac, h, two}];\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm, CellTags->"DoubleStep"], Cell[TextData[{ "Handler for Automatic values for the PrecisionGoal and AccuracyGoal \ options. This is essentially identical to how built-in ", StyleBox["Mathematica", FontSlant->"Italic"], " functions handle these values. " }], "Text"], Cell["\<\ AutoAGPG[wprec_, apgin_, pora_] := Module[{apg = apgin}, If[apg === Automatic, \t\tapg = If[wprec <= $MachinePrecision, 6, wprec - 10], \t\tIf[Not[NumberQ[apg] && Positive[apg]], \t\t\tIf[pora === Precision, \t\t\t\tMessage[SplitStepSolve::precg, apg], \t\t\t\tMessage[SplitStepSolve::accg, apg]]; \t\t\tReturn[$Failed]; \t\t] ]; SetPrecision[10.^(-apg), wprec] ];\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ The possibility of using an absorbing boundary potential is a property of the \ equation. The profile has been chosen to be smooth which helps minimize \ reflection. You need to pick an extent and size which balances reflection \ and transmission.\ \>", "Text"], Cell["\<\ AbsorbingBoundaryConditionPotential[g_, {xmin_, xmax_}, {xL_, xR_}][x_] := - I (g/2) If[x < xL, \t1 - Cos[Pi (x - xL)/(xmin - xL)], \tIf[x > xR, 1 - Cos[Pi (x - xR)/(xmax - xR)], 0]]\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ This is run after each set of steps. In addition to the values chosen, a rule \ Time\[Rule]t is always constructed for reference.\ \>", "Text"], Cell["\<\ StepOutput[t_, u_, things_, grid_] := Block[{datum, L, ua, plot = False}, \tL = n*(grid[[2]] - grid[[1]]); \tua = Abs[u]; \tdatum = Join[{Time->t},things]; \tDo[ \t\tdatum[[j]] = Switch[datum[[j]], \t\t\tValue, Value->u, \t\t\tAbsValue, AbsValue->Abs[u], \t\t\tMaxValue, MaxValue->Max[Abs[u]], \t\t\tMaxPosition, tmp = Abs[u]; MaxPosition->Position[tmp, Max[tmp]][[1,1]], \t\t\tGrid, Grid->grid, \t\t\tMass, Mass->Mass[ua, L], \t\t\tCentroid, Centroid->Centroid[ua, grid], \t\t\tVariance, Variance->Variance[ua, grid], \t\t\tMomentum, Momentum->Momentum[u, L, FourierDVector[1,Length[u], L]], \t\t\tPlot | {Plot,popts___}, \t\t\t\tpopts = Sequence[]; \t\t\t\tIf[ListQ[datum[[j]]] && Length[datum[[j]]] > 1, \t\t\t\t\tpopts = Apply[Sequence, Drop[datum[[j]],1]]]; \t\t\t\tplot = True;Show[Block[{$DisplayFunction = Identity}, { \t\t\t\t\tListPlot[Transpose[{grid, Abs[u]}], PlotJoined->True], \t\t\t\t\tListPlot[Transpose[{grid, Re[u]}], PlotJoined->True, \ PlotStyle->RGBColor[1,0,0]], \t\t\t\t\tListPlot[Transpose[{grid, Im[u]}], PlotJoined->True, \ PlotStyle->RGBColor[0,1,0]], \t\t\t\t\tListPlot[Transpose[{grid, Abs[u]}], PlotStyle->{PointSize[0.025], \ RGBColor[0,0,1]}]}], popts]; Plot \t\t], \t\t{j,2,Length[datum]}]; \tIf[TrueQ[plot],datum = DeleteCases[datum, Plot]]; \tdatum];\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm], Cell["\<\ Functions to compute output quantities, conserved quantities, etc.\ \>", "Text"], Cell["\<\ Mass[ua_, {xmin_, xmax_}] := Mass[ua, xmax-xmin] Mass[ua_, L_] := Tr[ua^2] L/Length[ua] Centroid[ua_, grid_] := Tr[ua^2 grid]/Tr[ua^2] Variance[ua_, grid_] := With[{c = Centroid[ua, grid]}, \tTr[(ua (c - grid))^2]/Tr[ua^2]] Momentum[u_, L_, fder_:Automatic] := Block[{fd = fder, ux}, \tIf[fder === Automatic, \t\tfd = FourierDVector[1,Length[u], L]]; \tux = InverseFourier[Fourier[u]*fd]; \t(L Im[Conjugate[u].ux]/Length[u]) ] End[]; (* Private *) EndPackage[];\ \>", "Input", PageWidth->Infinity, InitializationCell->True, ShowSpecialCharacters->False, FormatType->InputForm] }, Open ]] }, FrontEndVersion->"4.0 for Microsoft Windows", ScreenRectangle->{{0, 1024}, {0, 695}}, AutoGeneratedPackage->Automatic, WindowSize->{1016, 668}, WindowMargins->{{0, Automatic}, {Automatic, 0}}, StyleDefinitions -> "ArticleClassic.nb" ] (*********************************************************************** Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. ***********************************************************************) (*CellTagsOutline CellTagsIndex->{ "DoubleStep"->{ Cell[24531, 700, 773, 27, 391, "Input", InitializationCell->True, CellTags->"DoubleStep"]} } *) (*CellTagsIndex CellTagsIndex->{ {"DoubleStep", 29640, 869} } *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1739, 51, 41, 0, 72, "Title"], Cell[1783, 53, 89, 2, 27, "Text"], Cell[1875, 57, 313, 11, 121, "Input", InitializationCell->True], Cell[2191, 70, 154, 4, 31, "Input", InitializationCell->True], Cell[2348, 76, 979, 18, 95, "Text"], Cell[3330, 96, 3127, 74, 769, "Input", InitializationCell->True], Cell[6460, 172, 358, 6, 44, "Text"], Cell[CellGroupData[{ Cell[6843, 182, 145, 4, 31, "Input", InitializationCell->True], Cell[6991, 188, 54, 1, 30, "Output"] }, Open ]], Cell[7060, 192, 324, 7, 44, "Text"], Cell[7387, 201, 221, 7, 31, "Input", InitializationCell->True], Cell[7611, 210, 231, 4, 27, "Text"], Cell[7845, 216, 290, 10, 103, "Input", InitializationCell->True], Cell[8138, 228, 411, 8, 44, "Text"], Cell[8552, 238, 350, 12, 139, "Input", InitializationCell->True], Cell[8905, 252, 884, 14, 78, "Text"], Cell[9792, 268, 432, 8, 131, "Input"], Cell[10227, 278, 132, 3, 27, "Text"], Cell[10362, 283, 333, 12, 139, "Input", InitializationCell->True], Cell[10698, 297, 250, 4, 44, "Text"], Cell[10951, 303, 315, 11, 121, "Input", InitializationCell->True], Cell[11269, 316, 168, 3, 27, "Text"], Cell[11440, 321, 200, 6, 31, "Input", InitializationCell->True], Cell[11643, 329, 232, 4, 27, "Text"], Cell[11878, 335, 220, 7, 49, "Input", InitializationCell->True], Cell[12101, 344, 608, 11, 61, "Text"], Cell[12712, 357, 434, 15, 175, "Input", InitializationCell->True], Cell[13149, 374, 808, 15, 61, "Text"], Cell[13960, 391, 1276, 31, 481, "Input", InitializationCell->True], Cell[15239, 424, 106, 3, 27, "Text"], Cell[15348, 429, 251, 9, 85, "Input", InitializationCell->True], Cell[15602, 440, 720, 14, 61, "Text"], Cell[16325, 456, 6709, 209, 3487, "Input", InitializationCell->True], Cell[23037, 667, 1277, 25, 131, "Text"], Cell[24317, 694, 211, 4, 27, "Text"], Cell[24531, 700, 773, 27, 391, "Input", InitializationCell->True, CellTags->"DoubleStep"], Cell[25307, 729, 249, 6, 27, "Text"], Cell[25559, 737, 509, 18, 247, "Input", InitializationCell->True], Cell[26071, 757, 273, 5, 44, "Text"], Cell[26347, 764, 321, 9, 85, "Input", InitializationCell->True], Cell[26671, 775, 153, 3, 27, "Text"], Cell[26827, 780, 1422, 37, 535, "Input", InitializationCell->True], Cell[28252, 819, 90, 2, 27, "Text"], Cell[28345, 823, 600, 25, 373, "Input", InitializationCell->True] }, Open ]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)