数値的な例
Wolfram LibraryLink でアクセスできる関数を使うと，Wolfram言語の柔軟性と汎用性を維持しながら数値計算を最適化することが可能となる．既存の大規模な数値コードがある場合は，Wolfram Symbolic Transfer Protocol (WSTP)と LibraryLink の両方とも，Wolfram言語から駆動するコードを繋ぎ合わせるよい方法となる．一方数値計算を開発しているなら，Wolfram言語でプロトタイプを作成できるし，どこか一部が障害となっている場合はLibraryFunctionを使ってその部分の高速化を図ることができる．LibraryFunctionはWolfram言語の数値計算関数とも直接繋がっているため，別のサンプル点で評価されるコードを書くことで，計算が高速化できることもある．このチュートリアルでは，LibraryLink の後者2つの使用方法に焦点を当てる．
このチュートリアルで示される例題のソースはドキュメントパクレットにある．
source = FileNameJoin[{$InstallationDirectory, "SystemFiles", "Links", "LibraryLink", "LibraryResources", "Source"}];FilePrint[FileNameJoin[{source, "demo_numerical.c"}]]
基本的な二次関数の例
DLLEXPORT int parabola(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
mreal x, a, f;
if (Argc != 2) return LIBRARY_FUNCTION_ERROR;
x = MArgument_getReal(Args[0]);
a = MArgument_getReal(Args[1]);
f = x*x - a;
MArgument_setReal(Res, f);
return LIBRARY_NO_ERROR;
}
lf = LibraryFunctionLoad["demo_numerical", "parabola", {Real, Real}, Real]関数がロードされたら，プロットしたり，サンプルしたり等，通常の関数と同じように使える．
ContourPlot[lf[x, a], {x, -10, 10}, {a, 0, 100}]
Plot[x /. FindRoot[lf[x, a], {x, 1}], {a, 0, 2}]このような簡単な関数については，関数の評価のコストと比較して，Wolfram言語から結果を得たりWolfram言語に結果を得たりするオーバーヘッドが大きいため，LibraryFunctionを使うことについて実際の利点はない．
マンデルブロ(Mandelbrot)集合
一方，関数が評価するのにより高価な場合は，大きな利点があることがある．その一例がマンデルブロ集合の画像作成に使える関数である．
以下に示す関数mandelbrotは 
から始めて，列が有界でない
4の外に至るまで，
の繰返し数を計算するものである．
がmandelbrot_max_iterationsに等しくなるまで繰返しが続くと，点 
は集合内であると想定され，0が返される
static mint mandelbrot_max_iterations = 1000;
DLLEXPORT int mandelbrot(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
mint n = 0;
mcomplex z = {0.,0.};
mcomplex c = MArgument_getComplex(Args[0]);
mreal rr = mcreal(z)*mcreal(z);
mreal ii = mcimag(z)*mcimag(z);
while ((n < mandelbrot_max_iterations) && (rr + ii < 4)) {
mcimag(z) = 2.*mcreal(z)*mcimag(z) + mcimag(c);
mcreal(z) = rr - ii + mcreal(c);
rr = mcreal(z)*mcreal(z);
ii = mcimag(z)*mcimag(z);
iteration++;
}
if (n == mandelbrot_max_iterations)
n = 0;
MArgument_setInteger(Res, n);
return LIBRARY_NO_ERROR;
}
Cコードには複素演算が組み込まれていないため，数学演算は実数演算を使って実装されている．
mlf = LibraryFunctionLoad["demo_numerical", "mandelbrot", {Complex}, Integer]これで関数が複素平面上の点群を評価するWolfram言語の通常の関数のように使えるようになった．
Plot3D[mlf[x + I y], {x, -2., .5}, {y, -1.25, 1.25}, PlotPoints -> 100, ColorFunction -> "Rainbow"]関数は非常に不連続であるので，これは集合を可視化する最良の方法ではないだろう．より一般的で格段に速い方法は，通常の点の格子上で関数をサンプリングして画像を作成することである．
n = 501;Timing[samples = Table[mlf[x + I y], {y, -1.25, 1.25, 2.5 / (n - 1)}, {x, -2., .5, 2.5 / (n - 1)}]; ]Image[Unitize[samples], "Bit"]関数値から色チャンネルへのマップを作成しRasterを使うことでカラー画像が作れる．Imageを使うこともできるが，下に示すようにインタラクティブに使う場合はRasterの方が格段に速く表示できる．
colormap = Function[If[# == 0, {0., 0., 0.}, {# / 25., # / 25., 1.}]];
Graphics[Raster[Map[colormap, samples, {2}]], ImageSize -> 512]colormap = Function[If[# == 0, {0., 0., 0.}, Part[r, #]]] /. r -> RandomReal[1, {1000, 3}];
Graphics[Raster[Map[colormap, samples, {2}]], ImageSize -> 512]このようなもののゴールは，画像からできるだけ素早くデータを取り出すことであることがよくある．上の例では，C関数
が行う計算の主要部分は非常に素早く行われるが，残りの部分はWolfram言語コンパイラを使うと高速化できる．Wolfram言語コンパイラとLibraryFunctionは緊密に統合されており，これにより非常に効率的になる．
cmf = With[{f = mlf, c = colormap}, Compile[{{x, _Complex}}, c[f[x]], RuntimeAttributes -> Listable]];ここでWithを使うのは，Compileは自身の引数を評価せず，mlfとcolormapは使われるとすぐに外部評価によってランタイムで解決する必要のあるシンボルであると想定されるからである．それらを明示することで，コンパイラは可能な限り最も効率的な評価を行う．
Listable属性とはリストで与えられたいくつかの複素点について評価できるということである．
cmf[I Range[0., 1., .1]]これはつまり，色チャンネルデータを取得するのに必要なのは，この関数に必要な全サンプル点を与えるこどだけであるということである．これを素早く行う方法の一つは，コンパイル済み関数を書くということである．
cpoints = Compile[{{xmin, _Real}, {xmax, _Real}, {nx, _Integer}, {ymin, _Real}, {ymax, _Real}, {ny, _Integer}}, Table[x + I y, {y, ymin, ymax, (ymax - ymin) / (ny - 1)}, {x, xmin, xmax, (xmax - xmin) / (nx - 1)}]];mimage[{xmin_, xmax_, nx_}, {ymin_, ymax_, ny_}] :=
Raster[cmf[cpoints[xmin, xmax, nx, ymin, ymax, ny]]]AbsoluteTiming[Graphics[mimage[{-2., .5, n}, {-1.25, 1.25, n}]]]これらのコンパイル済み関数を使って画像全体を取得するのに必要な時間は，単にサンプル点で関数を評価するのより速くなる．これは上記サンプル点の評価に示される計算時間の多くはPackedArrayの代りにListとしてTableを構築するの使われているからである．コンパイル済み関数では全体を通してパックアレーが使われる．
複数のコアを持つマシンを使用しているなら，並列化によってさらに計算速度を増加させることができる．
cmf = With[{f = mlf, c = colormap}, Compile[{{x, _Complex}}, c[f[x]], RuntimeAttributes -> Listable, Parallelization -> True]];AbsoluteTiming[mimage[{-2., .5, n}, {-1.25, 1.25, n}];]高速化は顕著ではあるが，プロセッサ数に比例はしない．これはサンプル点の構築や画像の生成等，1つのプロセッサ上で行われる部分もあるからである．
このように速い更新では，集合のいくつかの側面が調べられるインタラクティブなインターフェースを設定しておくと便利である．EventHandlerおよびDynamicを使った設定はやや複雑であるが，効率的に動作させるようにする決め手は，上で設定したmimage関数のスピードである．
mandelPlot[{x0_, x1_, nx_}, {y0_, y1_, ny_}] := Panel[DynamicModule[{dragstart = {0, 0}, dragval = {0, 0}, box = False, getImage = True, getBB = False, img, pan = False, shiftBB = False, xlo, xhi, ylo, yhi, blo, bhi, toZ, shift, rescale},
rescale[{lo_, hi_}][s_] := Block[{c = (hi + lo) / 2, r = .5 10. ^ s}, {c - r, c + r}];
{xlo, xhi} = Sort[{x0, x1}];
{ylo, yhi} = Sort[{y0, y1}];
Column[{
Row[{"SubscriptBox[s, x]", Slider[Dynamic[RealExponent[xhi - xlo], Function[{xlo, xhi} = rescale[{xlo, xhi}][#];getBB = False; getImage = True]], {-10, 1}], Dynamic[RealExponent[xhi - xlo]]}],
Row[{"SubscriptBox[s, y]", Slider[Dynamic[RealExponent[yhi - ylo], Function[{ylo, yhi} = rescale[{ylo, yhi}][#];getBB = False; getImage = True]], {-10, 1}], Dynamic[RealExponent[yhi - ylo]]}], Button[Dynamic[If[TrueQ[pan], "Click to Zoom", "Click to Pan"]], pan = TrueQ[Not[pan]]], Dynamic[If[pan, "Drag in image to Pan", "Drag in image to select new frame"]],
EventHandler[Dynamic[
If[getImage,
If[getBB && Norm[dragval - dragstart] > 1.*^-10,
toZ = Function[{n, ll, ur}, With[{diff = ur - ll}, Function[ll + diff(# / n)]]][{nx, ny}, {xlo, ylo}, {xhi, yhi}];
If[pan,
shift = toZ[dragstart] - toZ[dragval];
{xlo, ylo} = blo + shift;
{xhi, yhi} = bhi + shift;
(* else *),
{xlo, ylo} = toZ[dragstart];
{xhi, yhi} = toZ[dragval];
{xlo, xhi} = Sort[{xlo, xhi}];
{ylo, yhi} = Sort[{ylo, yhi}];
];
getBB = False;
];
img = mimage[{xlo, xhi, nx}, {ylo, yhi, ny}];
getImage = False;
];
Show[Graphics[{img, If[TrueQ[box], {Opacity[.5], Yellow, Rectangle[dragstart, dragval]}, {}]}], ImageSize -> {nx, ny}]], "MouseDown" :> (dragval = dragstart = MousePosition["Graphics"]; If[pan, blo = {xlo, ylo};bhi = {xhi, yhi}, box = True]), "MouseDragged" :> (dragval = MousePosition["Graphics"]; getBB = pan; getImage = pan), "MouseUp" :> (box = False; getBB = True; getImage = True)]
}]]]マシンがかなり速ければ，各方向512ポイントでこれを実行することができる．そうでなくても200ポイントで適切な画像が得られる．子rを実行するためには
mandelPlot[{-2.5,.5,n},{-1.25,1.25,n}]を使うとよい．ここで n は各方向の点の数を示す．ズームの倍率をかなり大きくすると，2つの理由で鮮明さが失われる．その一つは，境界を鮮明にするにはCコードでエンコードされた1000回の繰返し上限以上が必要になるという点である．2つ目は場合によっては近隣の点の軌道を実際に解像するのに，より高度な精度が要求されることがあるという点である．
Wolfram言語コンパイラは，コンパイルするとLibraryFunctionとしてロードできる共有ライブラリになるCコードを自動的に生成する．例えば，生成されたCコードはより遅いが，手で書いた関数と比べるとそれほど遅くない．
mcf = Compile[{{c, _Complex}}, Module[{n = 0, z = 0. + 0. I, maxn = 1000},
While[n < maxn && Re[z] ^ 2 + Im[z] ^ 2 < 4, z = z ^ 2 + c;n++];If[n == maxn, 0, n]], RuntimeOptions -> "Speed", CompilationTarget -> "C"];cmcf = With[{f = Last[mcf], c = colormap}, Compile[{{x, _Complex}}, c[f[x]], RuntimeAttributes -> Listable, Parallelization -> True]];
mimagecf[{xmin_, xmax_, nx_}, {ymin_, ymax_, ny_}] :=
Raster[cmcf[cpoints[xmin, xmax, nx, ymin, ymax, ny]]]AbsoluteTiming[mimagecf[{-2., .5, n}, {-1.25, 1.25, n}];]ブリュセレータ(Brusselator)偏微分方程式
化学反応のモデル化のためのブリュセレータ偏微分方程式問題（[HNW93]と[HW96]に記載されている）は，記号前処理のコンポーネントをいくつか含むNDSolve等の関数で使われるLibraryFunctionを定義するよい例である．
α = 1 / 50;
PDE = {D[u[t, x], t] == 1 + u[t, x] ^ 2 v[t, x] - 4 u[t, x] + α D[u[t, x], {x, 2}],
D[v[t, x], t] == 3 u[t, x] - u[t, x] ^ 2 v[t, x] + α D[v[t, x], {x, 2}]};
BC = {u[t, 0] == u[t, 1] == 1, v[t, 0] == v[t, 1] == 3, u[0, x] == 1 + Sin[2π x], v[0, x] == 3};問題はNDSolveで直接解くことができる．この例の論点の一部は，NDSolveと部分的に手作業でコードを作成した離散化との解像の速度を比較するということである．
Timing[sol = NDSolve[{PDE, BC}, {u, v}, {t, 0, 10}, {x, 0, 1}]]Plot3D[u[t, x] /. sol, {t, 0, 10}, {x, 0, 1}]参考文献[HW96]では二次の差分と 
点の離散化を使って問題が設定されており，大きい 
については「硬い問題の全特徴を有している」と言っている．これと同じ離散化は，NDSolveを使って二次差分を強制し，線の方法のメソッドオプションを使って空間格子点の数を設定することで行える．NDSolve`ProcessEquationsを使い，結果のNDSolve`StateDataオブジェクトから関数を抽出してNDSolveで効率的に解き，常微分方程式の右辺を得ることもできる．
n = 200;state = First[NDSolve`ProcessEquations[{PDE, BC}, {u, v}, {t, 0, 10}, {x, 0, 1}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 2, "MaxPoints" -> n, "MinPoints" -> n}}]]init = First[state["SolutionVector"]];
rhs = state["NumericalFunction"];demo_numerical.cの関数brusselator_pde_rhsは同じ離散化が返されるように定義されている．
rhslf = LibraryFunctionLoad["demo_numerical", "brusselator_pde_rhs", {Real, {Real, 1}}, {Real, 1}]
は離散化された右辺に明示的には現れないが，関数はそれを第1引数に取るように定義されている．これによりNDSolveは非常に非効率的なインターフェースを作ることができる．この関数が与える値はNDSolveが生成した関数が与える値とほぼ同じである．
Max[Abs[rhs[0., init] - rhslf[0., init]]]わずかな差分は単に計算操作の順における小さな差のためである．
reps = 10 ^ 4;
{First[AbsoluteTiming[Do[rhs[0., init], {reps}]]], First[AbsoluteTiming[Do[rhslf[0., init], {reps}]]]}LibraryFunctionの方が約20倍速い．
NDSolveは偏微分方程式を解くのに，常微分方程式
の一次系の右辺に関数を使用した．ベクトルの初期条件を使うと，これはNDSolveに直接与えることができる．
AbsoluteTiming[vsol = NDSolve[{U'[t] == rhs[t, U[t]], U[0] == init}, U, {t, 0, 10}, AccuracyGoal -> 4, PrecisionGoal -> 4]]AbsoluteTiming[vsollf = NDSolve[{U'[t] == rhslf[t, U[t]], U[0] == init}, U, {t, 0, 10}, AccuracyGoal -> 4, PrecisionGoal -> 4]]LibraryFunctionは評価するのに20倍近く速くても，解には3倍もかかる．これはヤコビアン(Jacobian)が必要とする硬さと，ヤコビアンが疎であるということに起因する．この問題を回避する方法は2つある．最も簡単なのは，ヤコビアンの疎のパターンを与え，比較的少ない評価で有限の差分が使えるようにし，疎な線形代数が使えるようにすることである．
疎のパターンを得るのが難しいことはよくあり，また非零であろう要素をすべて得ることが重要である．以下では，右辺が 
に依存していないという事実を利用し，有限差分のみ使用する．
test = RandomReal[1, 2 n];
clf = rhslf[0., test];
jac = Transpose[SparseArray[Table[rhslf[0., test + .1 UnitVector[2 n, i]] - clf, {i, 2 n}]]];
jpat = jac["PatternArray"]AbsoluteTiming[vsol = NDSolve[{U'[t] == rhslf[t, U[t]], U[0] == init}, U, {t, 0, 10}, AccuracyGoal -> 4, PrecisionGoal -> 4, "Jacobian" -> {"FiniteDifference", "Sparse" -> jpat}]]解析的ヤコビアンを定義するのは，パターンを正しく得なくてはならない上に，そのパターンを使って値のベクトルを与える関数から疎配列を構築しなければならないので，1ステップ複雑になる．
jplf = LibraryFunctionLoad["demo_numerical", "brusselator_pde_jacobian_positions", {Real, {Real, 1}}, {Integer, 2}];jvlf = LibraryFunctionLoad["demo_numerical", "brusselator_pde_jacobian_values", {Real, {Real, 1}}, {Real, 1}];With[{pos = jplf[0., init]}, jfun[t_ ? NumberQ, U_] := SparseArray[pos -> jvlf[t, U]]];ここではjplfが一度だけ評価され，定義に位置が含まれるようにするためにWithが使われている．位置が時間やベクトルに依存する場合は，定義の中で明示的にjplf[t,U]を使えばよい．記号的評価を避けるためにPatternTest t_?NumberQが使われているが，これにより構成が正しくなくなる．
Max[Abs[jfun[0., init] - rhs["Jacobian"[0., init]]]]有限差分は機械精度の半分までしか正確でないので，ここは差の大きさが想定される．
AbsoluteTiming[vsol = NDSolve[{U'[t] == rhslf[t, U[t]], U[0] == init}, U, {t, 0, 10}, AccuracyGoal -> 4, PrecisionGoal -> 4, "Jacobian" -> jfun[t, U[t]]]]速度は有限差分を使った場合と同程度である．計算によってはヤコビアンの方が他のに比べてエラーになりやすいことがある．NDSolveについては，多くの場合ある程度近いヤコビアンで十分であることが多い．一方でFindRootやFindMinimumのような関数については，導関数の確度は解の確度に大きく関係する．
NDSolveを使って特定の値の 
に対する右辺関数のためのCコードを生成することもできる．メソッドオプションを"ExpandEquationsSymbolically" -> Trueとすると，ベクトル関数を使う代りに線の方法が系を記号的に拡張する．Compiledオプションでは，その中で与えられたオプションをコンパイラに渡す．
AbsoluteTiming[gstate = First[NDSolve`ProcessEquations[{PDE, BC}, {u, v}, {t, 0, 10}, {x, 0, 1}, Method -> {"MethodOfLines", "ExpandEquationsSymbolically" -> True, "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> 2, "MaxPoints" -> n, "MinPoints" -> n}}, Compiled -> {True, CompilationTarget -> "C"}]];
grhs = gstate["NumericalFunction"];]この処理を行うには長い時間がかかる．これは生成されたコードが非常に大きく，Cコンパイラがコードをコンパイルして最適化するのに時間がかかるからである．ヤコビアンが最初に評価されたときに，Cコードも生成されてコンパイルされるので，最初の評価には時間がかかる．
AbsoluteTiming[Max[Abs[jfun[0., init] - grhs["Jacobian"[0., init]]]]]AbsoluteTiming[vsol = NDSolve[{U'[t] == grhs[t, U[t]], U[0] == init}, U, {t, 0, 10}, AccuracyGoal -> 4, PrecisionGoal -> 4]]そのためコンパイルのための顕著な時間コストにおいて，生成されたコードの速度は手作業でコードした関数とほぼ同じである．生成されたコードはあまり柔軟ではないことに注意されたい．
や 
の別の値についてすべて再び生成しなければならない．
ダフィング(Duffing)方程式の折りたたみ
ダフィング方程式は定期的に駆動する固定梁であると考えることができる．これは非線形の復元力を持つ振動子を表す．
ϵ = 3 / 10; ω = 1;γ = 15 / 100;
deq = x''[t] + γ x'[ t] - x[t] + x[t] ^ 3 == ϵ Cos[ω t];これらのパラメータの値について，軌道は確定的にカオス的である．NDSolveが計算した解からこれを表示することができる．
pplot = ParametricPlot[Evaluate[{x[t], x'[t]} /. NDSolve[{deq, x'[0] == x[0] == 0}, x, {t, 0, 25}]], {t, 0, 25}]カオス軌道は一般に初期条件に非常に敏感である．この例では，初期条件の線について軌跡の変化をどのように追跡できるかが示されている．
これを行う一つの方法に，曲線を特定の時間までに解決できるように線に十分近い空間から始めるということである．この方法は比較的簡単だが計算的に集中し，さらにいくつの点で開始すべきかを予測するのが難しい．
別の方法は，軌跡が進むにつれて特徴が自動的に解決されるように適応的に点を追加するというものである．これはすべての軌跡を時間を増やしながら進め，ある基準に応じて点を追加して行う．この例で使われた基準は単に点と点の間の空間である．曲率等，より高度なテストも使えるが，単純に距離を使うと物事がより簡単にすむ．
解を前進させる最も効率のよい方法は，ベクトルの初期条件で記述された一次系でNDSolveを使うものである．
Clear[advance];advance[{initx_, inity_}, t0_, t1_, opts___] := Block[{t, x, y}, {x[t1], y[t1]} /. First[NDSolve[{x'[t] == y[t], x[t0] == initx, y'[t] == ϵ Cos[ω t] + x[t](1 - x[t]^2) - γ y[t], y[t0] == inity}, {x, y}, {t, t0, t1}, opts]]];
v = {0., 0.};
Δt = .1;
points = Table[v = advance[v, t - Δt, t], {t, Δt, 25, Δt}];
Show[pplot, Graphics[{Red, Point[points]}]]NDSolveでは初期条件をベクトルにすることができるため，advance関数は多くの軌跡を一度に扱うことができる．初期値 x と y は別々のベクトルなので，{{x1,y1},{x2,y2,…}を使う代りに{{x1,x2,…},{y1,y2,…}}を使って x と y の座標を分けたことに注意されたい．
v = {Range[0., 1., .1], ConstantArray[0., {11}]};
Δt = .1;
points = Table[Transpose[v = advance[v, t - Δt, t]], {t, Δt, 1.5, Δt}];
ListLinePlot[Transpose[points]]点と点の間の距離が 
において大きく増加したことが分かる．この微調整方法は距離が大きくなりすぎる前に新しい点を加えるというものである．
refine[dist_, {px_, py_}] := Module[{dx = Differences[px], dy = Differences[py], k, s, e, rx, ry},
s = MapThread[Function[
k = Max[1, Ceiling[Sqrt[(#1 ^ 2 + #2 ^ 2)] / dist]];
N[Range[0, k - 1]] / N[k]],
{dx, dy}];
rx = Apply[Join, Append[MapThread[#1 + #2 * #3&, {Drop[px, -1], dx, s}], px[[{-1}]]]];
ry = Apply[Join, Append[MapThread[#1 + #2 * #3&, {Drop[py, -1], dy, s}], py[[{-1}]]]];
{rx, ry}]advance関数にはいくつでも点が使えるため，指定の初期条件集合を前進させるためにしなければいけないことは，refineを使って点が十分に近いようにし，これらの点をすべて進め，それを繰り返すことである．
v = {{0., 1.}, {0., 0.}};
dist = 0.025;
Δt = 0.1;
points = Table[v = refine[dist, v];Transpose[v = advance[v, t - Δt, t]], {t, Δt, 1.5, Δt}];
ListPlot[points]
の値をより大きくすると，計算に長い時間がかかることがある．これは隣接する2点は指数関数的に発散することが多く，加えられる点の数も時間に対して指数関数的に増えるからである．
計算時間を減らすには，まずどこに多くの時間がかけられているかを見極める必要がある．この計算では，考えなければいけない要素は，微調整と解の前進の2つだけである．時間の多くがどこに使われているかを見るには，軌跡が延びるにしたがって各ステップが必要とする時間をプロットすればよい．
AbsoluteTiming[v = {{0., 1.}, {0., 0.}};
dist = 0.025;
Δt = 0.1;
Tend = 20;
ListLinePlot[Transpose[Table[{{t, Timing[v = refine[dist, v]][[1]]}, {t, Timing[v = advance[v, t - Δt, t]][[1]]}}, {t, Δt, Tend, Δt}]]]]プロットを見ると，各ステップの微調整部分で半分以上の時間が使われているのが分かる．したがって，その部分を早くすると，計算全体の速度に最も大きく影響があるであろうと考えられる．
微調整関数はパックアレーとしてコンパイルや表現ができない異なる長さの配列を扱うため，Wolfram言語ではスピードアップが難しい．したがって，この部分を行うためのLibraryFunctionを設定するのが理にかなっている．
refinelf = LibraryFunctionLoad["demo_numerical", "refine", {Real, {Real, 2}}, {Real, 2}];
AbsoluteTiming[v = {{0., 1.}, {0., 0.}};
dist = 0.025;
Δt = 0.1;
Tend = 20;
ListLinePlot[Transpose[Table[{{t, Timing[v = refinelf[dist, v]][[1]]}, {t, Timing[v = advance[v, t - Δt, t]][[1]]}}, {t, Δt, Tend, Δt}]]]]微調整にLibraryFunctionを使うと，計算が3倍早くなり，微分方程式の解の計算においては，計算時間が大半を占めることになる．
多くの場合，NDSolveに組み込まれたより高度な適応時間ステップ法を使って微分方程式を数値的に解くと，より効率的である．しかしこの例では，次の2つの事柄が組み合さって，簡単な固定ステップサイズのソルバを書くことが望ましい状態になっている．1つ目はメソッドによらず，局所誤差は最終的に指数関数的に大きくなるため，計算された解は実際の解を暗示しているにすぎず，長時間実行する場合は局所誤差を微調整しても，その誤差が大きすぎない限りは大きな違いは生じない．2つ目は，比較的高次の方法において，微調整間のステップは常微分方程式ソルバに対して適切であるということである．したがって，微調整はオーバーヘッドなしで非常に素早く結果を得るためのソルバの一ステップのみで置き換えることが可能である．ここに適したソルバは，古典的な4次のルンゲクッタ(Runge–Kutta)法であるが，より低い次数の方法で時間をさらに削減することも可能であるかもしれない．
advancelf = LibraryFunctionLoad["demo_numerical", "duffing_crk4", {{Real, _}, Real, Real}, {Real, _}]
AbsoluteTiming[v = {{0., 1.}, {0., 0.}};
dist = 0.025;
Δt = 0.1;
Tend = 20;
Do[v = advancelf[refinelf[dist, v], t - Δt, t], {t, Δt, Tend, Δt}]]単一ステップの前進LibraryFunctionは格段に速いことが分かる．時間の差は小さすぎてあまり意味を成さないので，プロットしなかった．さらに先の時間まで行き，時間の差をプロットするなら，時間の大半が微分方程式の解法に使われていることが分かるだろう．
繰返しが十分に最適化されたので，さらに大きな計算ができるようになった．比較的小さい列から始めて展開していくことで得られるプロットの列は，どのように近くの初期条件が伸ばされ折りたたまれていくかを示す．
AbsoluteTiming[v = {{0., 0.01}, {0., 0.}};
dist = 0.025;
Δt = 0.1;
ΔT = 1.0;
Tend = 60;
plots = Table[
Do[v = advancelf[refinelf[dist, v], t - Δt, t], {t, T - 1 + Δt, T, Δt}];
points = Transpose[v];
Graphics[Line[points, VertexColors -> Map[Hue, N[Range[Length[points]]] / Length[points]]], Axes -> True, PlotRange -> {{-1.5, 1.5}, {-1, 1}}],
{T, ΔT, Tend, ΔT }];]Take[plots, {10, 60, 10}]非常に小さな線分がやがて引き伸ばされ，折れ曲がって元に戻ってくることが分かる．これを見るのに適した方法は，ListAnimate[plots]を入力することである．サイズの都合上，出力はこのノートブックには含めない．