I wonder how can I implement dual numbers in Mathematica, so that all functions work well with them (as with complex numbers).
Particularly, for each function $f$, $f(\varepsilon)=f^\prime(0)\varepsilon$
Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this communityI wonder how can I implement dual numbers in Mathematica, so that all functions work well with them (as with complex numbers).
Particularly, for each function $f$, $f(\varepsilon)=f^\prime(0)\varepsilon$
Here's a tiny piece of advice I follow: any time I want to implement a new, exotic number system in Mathematica, the first thing I do is to look within the implementation of the Quaternions`
package, and try to adapt/emulate the constructions in the package to the number system I am trying to implement.
Having said this, here's a bunch of rules for doing basic arithmetic with dual numbers, as well as a general rule for evaluating functions with dual number arguments (notice the heavy use of TagSetDelayed[]
so that the rules are associated with Dual[]
and not the arithmetic/transcendental functions):
ScalarQ[c_] := ! MatchQ[Head[c], Dual] && NumericQ[c]
DualEpsilon = Dual[0, 1];
Dual /: Dual[a_, 0] := a;
Dual /: c_?ScalarQ + Dual[a_, b_] := Dual[c + a, b];
Dual /: Dual[a_, b_] + Dual[c_, d_] := Dual[a + c, b + d];
Dual /: c_?ScalarQ*Dual[a_, b_] := Dual[c a, c b];
Dual /: Dual[a_, b_]*Dual[c_, d_] := Dual[a c, b c + a d];
Dual /: Power[d_Dual, n_Integer?Positive] :=
Fold[(If[#2 == 1, d, 1] #1) #1 &, d, Rest[IntegerDigits[n, 2]]];
Dual /: Power[Dual[a_, b_], -1] := Dual[1/a, -b/a^2];
Dual /: Power[d_Dual, n_Integer?Negative] := 1/Power[d, -n];
Dual /: Power[Dual[a_, b_], x_?ScalarQ] := Dual[a^x, b x a^(x - 1)];
Dual /: Abs[Dual[a_, b_]] := Dual[Abs[a], b Sign[a]];
Dual /: Sign[Dual[a_, b_]] := Sign[a];
Dual /: f_[d__Dual] /; MemberQ[Attributes[f], NumericFunction] := With[{args = {d}},
Dual[f @@ args[[All, 1]],
(Derivative[##][f] @@ args[[All, 1]]) & @@@
IdentityMatrix[Length[args]].args[[All, 2]]]]
Samples:
Dual[a, b]/Dual[c, d]
Dual[a/c, b/c - (a d)/c^2]
Dual[a, b]^(2/3)
Dual[a^(2/3), (2 b)/(3 a^(1/3))]
Exp[Dual[a, b]]
Dual[E^a, b E^a]
3 Dual[a, b]^2 - 2 Dual[a, b]
Dual[-2 a + 3 a^2, -2 b + 6 a b]
Dual[a, b]^Dual[c, d]
Dual[a^c, a^(-1 + c) b c + a^c d Log[a]]
The last rule, while fairly general, will give erroneous results for functions like Floor[]
, Mod[]
, etc. Adding the rules needed for proper evaluation of those, as well as adding formatting rules (such that e.g. Dual[a, b]
prints as $a+b\varepsilon$ in StandardForm
) is left as an exercise.
Dual /: c_?ScalarQ * Dual[a_, b_] /; ScalarQ[c] := Dual[c a, c b];
?
$\endgroup$
Jun 11, 2015 at 21:52
If you use the matrix representation, addition and multiplication works by default, but you need to use matrix multiplication and not just space (which is element-wise multiplication), for instance f[a_,b_]:=a.b+b.b
would for instance work exactly like you would expect if a and b where dual numbers.
If you just want to have the rule enforce that you mention, you could just write that rule out:
dualE /: f_[dualE] := f'[0] dualE
You can add other rules to make your other algebra work out aswell for intance for multiplication:
dualE /: Power[dualE, 2] := 0
(a + dualE b) (c + dualE d) // Expand // Collect[#, dualE] &
Update
If you want the full function application rule, then there is a problem with just defining it, lets assume you wanted to write down: f_[a+b dualE]:=f[a]+b f'[0] dualE
if we expand this into the expression form it reads: f_[Plus[a,Times[b,dualE]]]:=...
now such a definition can't have it's pattern attached to dualE, since it doesn't appear in the top level. What you would then need is to create a symbol that represents a dual number with both parts included, as Sasha suggested in a comment. Here is an example writing up rules for Times
and Plus
:
dualE:=dualNumber[0,1];
Times[dualNumber[r1_,d1_],dualNumber[r2_,d2_]] ^:= dualNumber[r1 r2,d2 r1+d1 r2]
Times[dualNumber[r_,d_],n_] ^:= dualNumber[n r,n d]
Plus[dualNumber[r_,d_],n_] ^:= dualNumber[r+n,d]
Then you can write the definition for function application with respect to dualNumber
dualNumber /: f_[dualNumber[r_, d_]] := f[r] + d f'[r] dualE
And if you would like the numbers to be printed like 2+ 3 duelE
rather then as dualNumber[2,3]
. Then you can define a MakeBoxes
rule for it.
MakeBoxes[dualNumber[a_, d_], StandardForm] ^:=
RowBox[{MakeBoxes@a, "+", MakeBoxes@(Times[d, "dualE"])}]
Now then this works together to allow:
Sin[(4 + dualE) (2 + 3 dualE)]
Sin[8] + dualE 14
Sin[(a + b dualE) (c + d dualE)]
Sin[a c] + (b c + a d) dualE
Edited to correct error caught by Rahul Narain, also to improve incorrect output formating in symbolic expressions.
f[a+b dualE]
, you would have to implement it. Then however it gets problematic. Note my update to the answer.
$\endgroup$
dualNumber /: f_[dualNumber[r_, d_]] := f[r] + f'[r] d dualE
instead, to match the intuition of $\varepsilon$ being an infinitesimal.
$\endgroup$
Possibly a function that applies functions, by expanding to first order in the dualE
part.
dualfunc[func_, d_] := Normal[Series[func[d], {dualE, 0, 1}]]
dualfunc[Sin, 3 + dualE]
(* dualE Cos[3] + Sin[3] *)
This has the inconvenience of being awkward, since it would need to be used with basic arithmetic (but as noted in other responses, UpValues
on dualE
could perhaps address this). An advantage is that it should work in a reasonably consistent manner, and (I hope) not require much beyond what I did here; might need some special fault-handling if you encounter the likes of 1/dualE
.
Well,
$Post = MatrixFunction[
Function[ε, #], {{0, 0}, {1, 0}}] /. {{a_,
b_}, {0, a_}} ->
a + ε b /. {{a_, 0}, {b_, a_}} ->
a + ε b &;
After this one can use any formulas with dual unity ε
:
I^ε
Out= 1 + (I ε Pi)/2
Coefficient
, D
or Cases
to obtain the the second term that multiplies the epsilon.
$\endgroup$
Nov 13, 2022 at 22:52
You can map any dual number $a + b \epsilon$ into a $2 \times 2$ matrix: $$ a + b \epsilon \mapsto \begin{pmatrix} a & b \cr 0 & a \end{pmatrix} $$
DualNumber[a,b]
object and overload Times
, Plus
and Power
to work on them. If you need also elementary functions to work on it, then going the route with matrices might be less work.
$\endgroup$
This answer combines the best traits ftom my previous answer and the answer by Daniel Lichyblau.
The following line adds the symbol for dual unity the same way as if it was a pre-defined constant like imaginary unit:
$Post=Normal[Series[#, {ε, 0, 1}]]&;
One even can use it in equations.
UPDATE
Alternatively (and better):
$Post=(#/.ε->0)+ε(D[#, ε]/.ε->0)&