# Multicomplex Mathematics


```python
# Standard libraries
import sys

# Package from conda
import numpy as np
import sympy as sy

# The MultiComplex library
import pymcx
```

## Getting started

The function 
$$ f= x\sin(x) $$
has the derivative
$$ f' = x\cos(x)+\sin(x) $$

In order to evaluate this derivative with MultiComplex algebra, we define a function that takes multicomplex arguments, in this case the lambda function ``lambda z: z*np.sin(z)``, and we compare the derivative obtained with the exact solution.

The function ``diff_mcx1`` takes as first argument the callable function of a single variable to be differentiated numerically, followed by the real value at which the derivative(s) are to be taken, and finally the number of derivatives that are to be taken.  All derivatives from order 1 to nderiv are returned in a list


```python
print(pymcx.diff_mcx1.__doc__)
```

    diff_mcx1(*args, **kwargs)
    Overloaded function.
    
    1. diff_mcx1(f: Callable[[MultiComplex<double>], MultiComplex<double>], x: float, numderiv: int, and_val: bool = False) -> List[float]
    
    2. diff_mcx1(f: Callable[[MultiComplex<double>], Tuple[MultiComplex<double>, MultiComplex<double>]], x: float, numderiv: int, and_val: bool = False) -> Tuple[List[float], List[float]]
    



```python
x = 0.1234
exact = x*np.cos(x) + np.sin(x)
nderiv = 4
and_val = False
mcx = pymcx.diff_mcx1(lambda z: z*np.sin(z), x, nderiv, and_val)[0]
error = mcx - exact
error
```




    0.0



## Derivatives of a function of one variable

Here we test some more interesting functions, and we show that the relative errors of these functions are all small, close to the numerical precision


```python
def deriver(var, func, nderiv):
    # The python translation of the sympy function
    # that takes numpy arguments (which allows for direct 
    # use of pymcx arguments)
    py_func = sy.lambdify(var, func, 'numpy')

    # Calculate "exact" derivatives from sympy and mpmath
    exacts = []
    for N in range(1, nderiv+1):
        f = sy.lambdify(var,sy.diff(func,var,N),'mpmath')
        exacts.append(f(0.1234))
    exacts = np.array(exacts)

    # Calculate MultiComplex derivatives for derivatives of order 
    # 1 to nderiv (inclusive) in one shot
    mcx = pymcx.diff_mcx1(py_func, 0.1234, nderiv)
    error = np.array(mcx - exacts)/np.array(exacts)
    return error

# Calculate some derivatives of functions to test out our implementation
y = sy.symbols('y')
display(deriver(y, 1.0/sy.log(y), 3))
display(deriver(y, sy.cos(y)*sy.sin(y)*sy.exp(y), 3))
display(deriver(y, 1/(4+sy.cos(y)*sy.sin(y)*sy.exp(y)/sy.cosh(y))-sy.log(y), 3))
```


    array([mpf('2.3990722300687498e-16'), mpf('5.7028372834229652e-15'),
           mpf('-2.440695170029012e-15')], dtype=object)



    array([mpf('1.7975675189132607e-16'), mpf('1.247770572193306e-16'),
           mpf('0.0')], dtype=object)



    array([mpf('2.1730817635703391e-16'), mpf('2.1651710041894285e-16'),
           mpf('-2.1372587912639966e-16')], dtype=object)


## Derivatives of multivariate function

In the multivariate case, the analysis is very similar.  We define a function that takes a _vector_ of multicomplex arguments, and in the diff_mcxN function specify the real values at which the derivatives are to be taken, and the number of times to take derivatives with respect to each of the independent variables


```python
def func(zs):
    x, y, z = zs
    return np.cos(x)*np.sin(y)*np.exp(z)
xs = [0.1234, 20.1234, -4.1234]
orders = [1, 1, 2]
def exact(zs):
    x, y, z=zs
    return -np.sin(x)*np.cos(y)*np.exp(z)

pymcx.diff_mcxN(func, xs, orders), exact(xs)
```




    (-0.0005830791997394103, -0.0005830791997394103)




```python
# A bad set of inputs (order too short), and an error returned
pymcx.diff_mcxN(func, xs, [1])
```


```python
def func(zs):
    x, y, z = zs
    return np.cos(x)*np.sin(y)*np.exp(z)

def der_sympy(xs):
    """ "Exact" solution in higher precision mathematics with sympy """
    x, y, z = sy.symbols('x,y,z')
    func = sy.cos(x)*sy.sin(y)*sy.exp(z)
    f = sy.diff(sy.diff(sy.diff(func,x),z,2),y)
    func = sy.lambdify((x,y,z),f,'mpmath')
    return func(*xs)

zs = []
h = 1e-50
L = 2
xs = [0.1234, 20.1234, -4.1234]
orders = [1, 1, 2]
numderiv = sum(orders)

# This block here is the implementation of the code that goes into
# diff_mcxN, and a simplified variant of it goes into diff_mcx1
for i in range(3):
    c = np.zeros((2**numderiv,))
    c[0] = xs[i]
    assert(orders==[1,1,2]) # The orders are hard-coded here..
    # <magic>
    if i == 0:
        c[1] = h
    elif i == 1:
        c[2] = h
    elif i == 2:
        c[4] = h
        c[8] = h
    # </magic>
    zs.append(pymcx.MultiComplex(c))

# Three values are the derivatives from multicomplex, derivative from symbolic math, 
# and derivative from N-dimensional derivative function
func(zs).get_coef()[2**L-1]/h**L, float(der_sympy(xs)), pymcx.diff_mcxN(func, xs, orders)
```




    (-0.0005830791997394104, -0.0005830791997394103, -0.0005830791997394103)



## Time Profiling


```python
x = pymcx.MultiComplex([1,2])
%timeit x*x
x = pymcx.MultiComplex([1,2,3,4])
%timeit x*x
nderiv = 4
%timeit pymcx.diff_mcx1(lambda x: x*np.sin(x), 0.1234, nderiv)
%timeit pymcx.diff_mcx1(lambda x: x*np.sin(x)*np.cos(x), 0.1234, nderiv)
```

    1.34 µs ± 6.95 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
    1.43 µs ± 10.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
    18 µs ± 55.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
    30.8 µs ± 90.7 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)


808a2faa1bdb3ae13d9bd0bf5b0270a251e1d0ef
ipynb
Jupyter Notebook
externals/mcx/Tutorial.ipynb
usnistgov/potter
8d8e396fb37caee6d9be37396c5a3b783d7ed0b1
MIT
2020-05-09T02:36:27.000Z
2021-01-12T00:46:54.000Z
externals/mcx/Tutorial.ipynb
svenpohl1991/potter
6c006561d0adc07d9458e567ef3eb3345de63340
MIT
externals/mcx/Tutorial.ipynb
svenpohl1991/potter
6c006561d0adc07d9458e567ef3eb3345de63340
MIT
2020-03-24T17:17:13.000Z
2022-03-03T17:43:39.000Z
Qwen/Qwen-72B
1. YES
2. YES