```python
import numpy as np
import scipy as sp
import sympy as smp
import matplotlib.pyplot as plt
from scipy.misc import derivative
```

# Symbolic Case

You know the formula explicitly but don't really feel like using a pencil to take the derivative

$$f(x) = ... $$

**Example:** $$f(x) = e^{-a\sin(x^2)}\cdot \sin(b^x) \cdot \ln(c\sin^2(x)/x)$$

First Define all your symbols in sympy


```python
x, a, b, c = smp.symbols('x a b c', real=True)
```


```python
x**2+smp.exp(a)
```




$\displaystyle x^{2} + e^{a}$




```python
x, a, b, c = smp.symbols('x a b c', real=True)
f = smp.exp(-a*smp.sin(x**2)) * smp.sin(b**x) * smp.log(c*smp.sin(x)**2 /x)
```

Compute derivatives using `smp.diff(f, x)` where $f$ is the function you want to take the derivative of and $x$ is the variable you are taking the derivative with respect to.


```python
dfdx = smp.diff(f, x)
dfdx
```




$\displaystyle - 2 a x e^{- a \sin{\left(x^{2} \right)}} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \sin{\left(b^{x} \right)} \cos{\left(x^{2} \right)} + b^{x} e^{- a \sin{\left(x^{2} \right)}} \log{\left(b \right)} \log{\left(\frac{c \sin^{2}{\left(x \right)}}{x} \right)} \cos{\left(b^{x} \right)} + \frac{x \left(\frac{2 c \sin{\left(x \right)} \cos{\left(x \right)}}{x} - \frac{c \sin^{2}{\left(x \right)}}{x^{2}}\right) e^{- a \sin{\left(x^{2} \right)}} \sin{\left(b^{x} \right)}}{c \sin^{2}{\left(x \right)}}$



Can take the nth derivative $d^n f/dx^n$ by putting the optional argument at the end `smp.diff(f,x,n)`.


```python
d4fdx4 = smp.diff(f, x, 4)
#d4fdx4
```

Can then compute numerical values by plugging in numbers


```python
d4fdx4.subs([(x,4),(a,1),(b,2),(c,3)]).evalf()
```




$\displaystyle -73035.8044625845$



Can also convert to a numerical function for plotting


```python
d4fdx4_f = smp.lambdify((x,a,b,c), d4fdx4)
```

Define $x$ and $y$ arrays using the numerical function above


```python
x = np.linspace(1,2,100)
```


```python
x = np.linspace(1,2,100)
y = d4fdx4_f(x, a=1, b=2, c=3)
```

Plot


```python
plt.plot(x,y)
plt.ylabel('$d^4 f / dx^4$', fontsize=24)
plt.xlabel('$x$', fontsize=24)
```

# Numerical Case

Your given data like such


```python
x, y = np.loadtxt('../data/sample_data1.txt')
```


```python
plt.plot(x, y, 'o--')
```

There are a few ways to go about this.

**1: The basic way**


```python
dydx = np.gradient(y,x)
```


```python
plt.plot(x,y, 'o--', label='$y(x)$')
plt.plot(x,dydx, 'o--', label='$y\'(x)$')
plt.legend()
```

This way works fine if the data is smooth but not if the data is noisy:


```python
x, y = np.loadtxt('../data/sample_data2.txt')
dydx = np.gradient(y,x)
```

Plot the noisey data


```python
fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
[a.legend() for a in ax]
plt.show()
```

Noise gets amplified in the derivative! This leads us to technique 2:

**2. The Clever Way**

Smooth data then take derivative. Consider the following new covid cases per day data


```python
x, y = np.loadtxt('../data/coviddata.txt')
dydx = np.gradient(y,x)
```

Taking the derivative naively gives a bad result


```python
fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
[a.legend() for a in ax]
plt.show()
```

Smooth the data by convolving it with a rectangle


```python
filt = np.ones(15)/15
```


```python
y_smooth = np.convolve(y, filt, mode='valid')
dysdx = np.gradient(y_smooth, x[7:-7])
```

Plot


```python
fig, ax = plt.subplots(1, 2, figsize=(10,3))
ax[0].plot(x,y, label='$y(x)$')
ax[0].plot(x[7:-7], y_smooth, label=r'$y_{{smooth}}(x)$')
ax[1].plot(x,dydx, label='$y\'(x)$', color='r')
ax[1].plot(x[7:-7],dysdx, label='$y_{smooth}\'(x)$', color='purple')
ax[1].set_ylim(-100,120)
ax[1].grid()
[a.legend() for a in ax]
[a.set_xlabel('Time [Days]') for a in ax]
ax[0].set_ylabel('Cases per Day')
ax[1].set_ylabel('$\Delta$ (Cases per Day) / $\Delta t$')
fig.tight_layout()
plt.show()
```

# Quasi-Symbolic Case

In this case you know your function $f(x) = ...$ but the function is not given by a typical expression. For example.

$$f(u) = \text{max} \left\{ \left|e^{-x_iu^2} - y_i\right| \right\}  $$

for some array of $x_i$'s and $y_i$'s


```python
x = np.linspace(0, 1, 500)
y = np.exp(-x*2.15**2) + 0.1*np.random.randn(len(x))
```


```python
plt.scatter(x,y)
plt.xlabel('$x_i$', fontsize=20)
plt.ylabel('$y_i$', fontsize=20)
plt.show()
```

Define function


```python
def f(u):
    return max(np.abs(np.exp(-x*u**2) -y))
```

Compute the function values


```python
u = np.linspace(0,10,40)
f_u = np.vectorize(f)(u)
```

Plot


```python
plt.plot(u, f_u, 'o--')
plt.xlabel('$u$', fontsize=20)
plt.ylabel('$f(u)$', fontsize=20)
plt.show()
```

You could always just take the derivative of the numerical array `f_u`, but there is a specific derivative function better assigned for this

Compute derivative


```python
dfdu = np.vectorize(derivative)(f, u, dx=1e-6)
```

Plot


```python
plt.plot(u, dfdu)
```

e3422be716c5f54dc015a37bfb46ba80b187d14c
ipynb
Jupyter Notebook
Maths/math3.ipynb
Sumanshekhar17/Pseudo-Spectral-Method
5d1c2449bafec89006de95f785fe569d488792d9
MIT
2021-09-05T18:21:22.000Z
2021-09-05T18:26:21.000Z
Maths/math3.ipynb
Sumanshekhar17/Pseudo-Spectral-Method
5d1c2449bafec89006de95f785fe569d488792d9
MIT
2021-09-05T18:26:10.000Z
2021-09-14T08:54:49.000Z
Maths/math3.ipynb
Sumanshekhar17/Pseudo-Spectral-Method
5d1c2449bafec89006de95f785fe569d488792d9
MIT
2021-09-06T16:36:35.000Z
2021-09-06T16:36:35.000Z
Qwen/Qwen-72B
1. YES
2. YES