Documentation Center

  • Trial Software
  • Product Updates

Integrate Custom Functions into MuPAD

MuPAD® provides a variety of tools for handling built-in mathematical functions such as sin, cos, exp, and so on. These tools implement the mathematical properties of these functions. Typical examples are the float conversion routine, the diff differentiation function, or the expand function, which you use to manipulate expressions:

float(sin(1));
diff(sin(x), x, x, x);
expand(sin(x + 1))

You can say that the mathematical knowledge about the built-in functions is distributed over several system functions: float knows how to compute numerical approximations of the sine function, diff knows the derivative of the sine function, and expand knows the addition theorems of the trigonometric functions.

When you implement your own function, you can integrate it into the MuPAD system, making other functions aware how to work with this new function correctly. If the new function consists only of built-in MuPAD functions, then you do not need to take extra steps. All MuPAD functions interact correctly with the new function:

f := x -> (x*sin(x)):
diff(f(x), x)

However, if you implement a function that is not composed of the standard MuPAD objects (for example, a new special function), you must distribute the knowledge about the mathematical meaning of the new function to standard MuPAD functions, such as diff, expand, float, and so on. This extra task is necessary for integrating the new function with the rest of the system. For example, you might want to differentiate an expression that contains both the new function and some built-in functions, and such differentiation is only possible via the MuPAD differentiation routine. Therefore, this routine must know how to handle the new symbol.

MuPAD uses function environments (domain type DOM_FUNC_ENV) to integrate functions into the system. A function environment stores special function attributes (slots) in an internal table. Whenever an overloadable system function, such as diff, expand, or float, encounters an object of type DOM_FUNC_ENV, it searches the function environment for a corresponding slot. If a system function finds the appropriate slot, it calls that slot and returns the value produced by the slot. All built-in MuPAD functions are implemented as function environments:

domtype(sin), domtype(exp)

You can call a function environment as you would call any MuPAD function or procedure:

sin(1.7), exp(1.0)

Suppose you implement the complete elliptic integral functions of the first and second kind, K(z) and E(z). These functions appear in different contexts, such as calculating the perimeter of an ellipsis, the gravitational or electrostatic potential of a uniform ring, and the probability that a random walk in three dimensions ever goes through the origin. The elliptic integrals have the following special values:

, E(1) = 1, , .

MuPAD provides the built-in functions ellipticE and ellipticK for computing these elliptic integrals. However, you can implement your own functions for the same task. For example, write the procedures ellipE and ellipK. These procedures define the values of the elliptic integrals for special values of x. For all other argument values, the values of elliptic integrals are unknown, and the procedures return the symbolic expressions ellipE(x) and ellipK(x). Use procname to return symbolic expressions:

ellipE :=
proc(x) begin
  if x = 0 then PI/2
  elif x = 1 then 1
  else procname(x) end_if
end_proc:
ellipK :=
proc(x) begin
  if x = 0 then PI/2
  elif x = 1/2 then 8*PI^(3/2)/gamma(-1/4)^2
  elif x = -1 then gamma(1/4)^2/4/sqrt(2*PI)
  else procname(x) end_if
end_proc:

ellipE and ellipK return special values for particular arguments. For all other arguments, they return symbolic expressions:

ellipE(0), ellipE(1/2),
ellipK(12/17), ellipK(x^2 + 1)

The first derivatives of these elliptic integrals are as follows:

, .

The standard MuPAD differentiation function diff does not know about these rules. Therefore, trying to differentiate ellipE and ellipK simply returns the symbolic notations of the derivatives:

diff(ellipE(x), x),
diff(ellipK(x), x)

To make diff work with the new functions, create function environments from the procedures ellipE and ellipK. In addition, function environments let you control the appearance of the symbolic function calls in outputs.

A function environment consists of three operands.

  • The first operand is a procedure that computes the return value of a function call.

  • The second operand is a procedure for printing a symbolic function call on the screen.

  • The third operand is a table that specifies how the system functions handle symbolic function calls.

To create function environments, use funcenv. For example, create function environments ellipE and ellipK. Use the second argument to specify that symbolic calls to ellipE and ellipK must appear as E and K outputs:

output_E := f -> hold(E)(op(f)):
ellipE := funcenv(ellipE, output_E):
output_K := f -> hold(K)(op(f)):
ellipK := funcenv(ellipK, output_K):

Although ellipE and ellipK are now function environments, you can call them as you would call any other MuPAD function:

ellipE(0), ellipE(1/2),
ellipK(12/17), ellipK(x^2+1)

The third argument funcenv is a table of function attributes. It tells the system functions (such as float, diff, expand, and so on) how to handle symbolic calls of the form ellipE(x) and ellipK(x). You can update this table specifying the rules for the new function. For example, specify the new differentiation rules by assigning the appropriate procedures to the diff slot of the function environments:

ellipE::diff :=
proc(f,x)
local z;
begin
z := op(f);
(ellipE(z) - ellipK(z))/(2*z) * diff(z, x)
end_proc:
ellipK::diff :=
proc(f,x)
local z;
begin
z := op(f);
(ellipE(z) - (1-z)*ellipK(z))/
(2*(1-z)*z) * diff(z, x)
end_proc:

Now, whenever f = ellipE(z), and z depends on x, the call diff(f, x) uses the procedure assigned to ellipE::diff:

diff(ellipE(z), z);
diff(ellipE(y(x)), x);
diff(ellipE(x*sin(x)), x)

The new differentiation routine also finds higher-order derivatives:

diff(ellipE(x), x, x)

Since the taylor function internally calls diff, the new differentiation routine also lets you compute Taylor expansions of the elliptic integrals:

taylor(ellipK(x), x = 0, 6)

If a derivative of a function contains the function itself, the integration routine has a good chance of finding symbolic integrals after you implement the diff attributes. For example, int now computes the following integrals:

int(ellipE(x), x)

int(ellipK(x), x)

Was this topic helpful?