how to do this kind of substitution in mathematica

Question

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.

Answer 1

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.

Answer 2

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 want g to be rewritten into something else
  before f "sees" it, but only if f does not hold its argument in position where g[something] stands

• Add an UpValue for g, with respect to expression f[g[_]] or similar. This will work if you want to alter the behavior of f in a soft way, that is, only on arguments of form matching the UpValue definition (f[g[_]] in this case). For this to work, either g should not have DownValues matching g[something] (so that it does not evaluate), or f should hold its argument. This is a typical use case when f is system function and g is user-defined, since this is the softest way to alter the behavior of the system function - note that while f behaves differently, the rule is attached to g.

• Add a DownValue to f. This is a "hard" way. Will work either if f holds its argument, or if g[something] will not by itself evaluate to something else, or you will have to temporarily add Hold - attribute to f, which is even more intrusive. If f is a system function, then this method is a last resort. Especially avoid this if f 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.