how to do this kind of substitution in mathematica
In case that I don't want to use mma's Erfc or Erf functions, I want to define my 开发者_如何学运维own
myErf[x_] := Integrate[Exp[-y^2/2], {y, -Infinity, x}];
And then I want mma to substitute any evalulation of integral with myErf[x], for example:
Integrate[Exp[-y^2+y/2], {y, x, Infinity}]
How to do this? Or in general, how to ask mma to evaluate an expression using the function I'd like it to use, instead of its own built in one? Is there any way of doing this?
Many thanks.
What I meant in my comment is that the naive
Unprotect[Integrate];
Integrate[Exp[-y_^2], {y_, -\[Infinity], x_}] := myErf[x]
to try to intercept some Integrate
evaluations does not work.
Actually it does work
In[]:= Integrate[Exp[-y^2],{y,-Infinity,z}]
Out[]= myErf[z]
(that will teach me for answering a SO question at midnight).
So maybe you don't need the myIntegrate
function given below.
Anyway... The two alternatives given in the comment are:
1) You could write your own Erf
integrator using something like
In[1]:= myErf[x_?NumericQ]:=NIntegrate[Exp[-y^2],{y,-\[Infinity],x}]
(* sames as Erf up to constant and scaling *)
In[2]:= Clear[myIntegrate]
myIntegrate[a_. Exp[b_. y_^2], {y_, -\[Infinity], x_}] /; FreeQ[{a, b}, y] := ConditionalExpression[myErf[x], b < 0]
myIntegrate[a_. Exp[b_. y_^2], {y_, xx_, x_}] := myIntegrate[a Exp[b y^2], {y, -\[Infinity], x}] - myIntegrate[a Exp[b y^2], {y, -\[Infinity], xx}]
myIntegrate[a_ + b_, lims__] := myIntegrate[a, lims] + myIntegrate[b, lims] (* assuming this is ok *)
myIntegrate[expr_, lims__] := Integrate[expr, lims]
In[8]:= myIntegrate[a Exp[x]+b Exp[-b x^2],{x,2y,z}]
Out[8]= ConditionalExpression[a (-E^(2 y)+E^z) - myErf[2 y] + myErf[z], b > 0]
Where at the end, the remaining integrals are passed to Integrate
. But you'd have to be smarter in the pattern matching (maybe including some transformations in intermediate steps) in order to catch all Erf
-like functions.
2)
You overwrite the Erf
-type functions so that they never evaluate using the built-in routines:
In[9]:= Clear[myErf];
Unprotect[Erf,Erfc,Erfi];
In[10]:= Erf[x__] := myErf[x]
Erfc[x__]:= myErfc[x]
Erfi[x__]:= myErfi[x]
In[13]:= {Erf[1], Erfc[2], Erfi[\[Pi]], Integrate[E^-x^2,x]}
Out[13]= {myErf[1], myErfc[2], myErfi[\[Pi]], 1/2 Sqrt[Pi] myErf[x]}
Where, obviously, you choose your own conventions for the myErf
functions.
Edit:
The simplest and maybe safest/softest method (thanks to Leonid) of doing what you want is to just Block[{Erf=...},...]
and let the Mathematica integrator do all of the work.
For example you myErf
was given as
In[1]:= myErf[x] == Integrate[Exp[-y^2/2], {y, -Infinity, x}]
Out[1]= myErf[x]==Sqrt[Pi/2] (1+Erf[x/Sqrt[2]])
So to solve the integral in your original post,
In[2]:= Block[{Erf = (Sqrt[2/Pi]*myErf[Sqrt[2] #] - 1 &)},
Integrate[Exp[-y^2 + y/2], {y, x, Infinity}]]
Out[2]= (E^(1/16)(Sqrt[\[Pi]] (-1 + 4 x + Sqrt[(-1 + 4 x)^2]) +
Sqrt[2] (1 - 4 x) myErf[Sqrt[2] Sqrt[(-(1/4) + x)^2]])
)/(2 Sqrt[(1 - 4 x)^2])
Edit 2:
Maybe you just need to be able to translate to and from myErf
expressions.
This means that it won't be automatic, but it will be flexible. Try
In[1]:= toMyErf=(FunctionExpand[#]/.Erf:>(Sqrt[2/Pi]*myErf[Sqrt[2] #]-1&)&);
fromMyErf=(#/.myErf:>(Sqrt[Pi/2] (1+Erf[#/Sqrt[2]])&)&);
In[3]:= Integrate[Exp[-y^2+y/2],{y,x,Infinity}]//toMyErf
Out[3]= 1/2 E^(1/16) Sqrt[\[Pi]] (2-Sqrt[2/\[Pi]] myErf[Sqrt[2] (-(1/4)+x)])
In[4]:= D[x*myErf[x+x^2],x]//fromMyErf//toMyErf
Out[4]= E^(-(1/2) (x+x^2)^2) x (1+2 x)+myErf[x+x^2]
Note that the FunctionExpand
command is there to rewrite Erfc[x]
as 1-Erf[x]
. You might want to implement that part slightly better.
A softer way is to use UpValues
:
ClearAll[withMyErf];
SetAttributes[withMyErf, HoldAll];
withMyErf[code_] :=
Block[{Exp},
Exp /: Integrate[Exp[(a : _ : 1)*x_^2], {x_, t0_, t_}] :=
ConditionalExpression[(Sqrt[Pi]*(myErf[Sqrt[-a]*t] - myErf[Sqrt[-a]*t0]))/(2*Sqrt[-a]),
Re[a] < 0];
code]
Example:
In[4]:= withMyErf[Integrate[Exp[-3*x^2], {x, -Infinity, t}]]
Out[4]= 1/2 Sqrt[\[Pi]/3] (myErf[Sqrt[3] t] - myErf[-\[Infinity]])
outside of our custom "dynamic scoping construct", everything will work as usual. Put the code where you want this to happen inside it. And, no changes ever made to Integrate
, which can be dangerous. All we risk here is that Exp will not evaluate, but most symbolic transformation rules for it are anyway local rules, I guess.
Edit:
To answer the general question, if you want to alter the standard way Mathematica evaluates things, you have several choices. Some of them have been mentioned in the previous answers, I just want to summarize them on a little more abstract level. Given an expression f[g[something]]
, typically you can:
Redefine DownValues for
g
- this will work if you wantg
to be rewritten into something else
beforef
"sees" it, but only iff
does not hold its argument in position whereg[something]
standsAdd an
UpValue
forg
, with respect to expressionf[g[_]]
or similar. This will work if you want to alter the behavior off
in a soft way, that is, only on arguments of form matching theUpValue
definition (f[g[_]]
in this case). For this to work, eitherg
should not haveDownValues
matchingg[something]
(so that it does not evaluate), orf
should hold its argument. This is a typical use case whenf
is system function andg
is user-defined, since this is the softest way to alter the behavior of the system function - note that whilef
behaves differently, the rule is attached tog
.Add a
DownValue
tof
. This is a "hard" way. Will work either iff
holds its argument, or if g[something] will not by itself evaluate to something else, or you will have to temporarily add Hold - attribute tof
, which is even more intrusive. Iff
is a system function, then this method is a last resort. Especially avoid this iff
is a non-trivial or frequently used and important built-in function, because it is hard to anticipate all the consequences of such modifications, while many system functions are using other system functions on the top level.Write your own lexical scoping constructs, which will analyze the code in unevaluated form and do the proper substitutions, and only then allow the code to evaluate. This will require that you have access to your code in one piece (that is, the code must exist in the form you want from the start), which is a severe restriction, so in practice this option is ruled out in most cases.
There is also a Block
trick, which allows you to "localize the damage" by making such redefinitions dynamically local to a block of code where you want them to apply. In some cases, I find it indispensable, because it allows ways to non-locally control the evaluation process which are not easily reproduced by other language constructs.
精彩评论