```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy
from scipy.integrate import odeint
from scipy.integrate import quad
from scipy.optimize import minimize
from scipy.optimize import fsolve
from scipy.optimize import least_squares
from scipy.interpolate import interp1d
```

Go slow, write the plan down, and check energy


```python
mpl.rc('figure',dpi=250)
mpl.rc('text',usetex=True)
```


```python
def add_labels(xlabel, ylabel, title):
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)
    plt.legend()
```


```python
def generate_samples(func, bounds, N): 
    """
    Assignment 5: function for sampling from a custom 1d distribution.
    
    Input:  func - any function
            bounds - tuple (a, b)
    
    Output: 1d samples of size N
    """
    a, b = bounds
    x = np.linspace(a, b)
    
    f_max = max(func(x))
    
    x_rand = np.random.uniform(a, b, size=N)
    y_rand = np.random.uniform(0, f_max, size=N)
        
    samples = x_rand[y_rand < func(x_rand)]

    return samples
```


```python
def spherical_to_xyz(r, theta, phi):
    "Generate a vector in xyz coordinates from r, theta, phi"
    
    return np.array([r*np.sin(theta)*np.cos(phi), r*np.sin(theta)*np.sin(phi), r*np.cos(theta)])
```

Solves a problem of the form

\begin{equation}
\frac{dy}{dt} = f(y, t)
\end{equation}



```python
def main(y, t):
    "Main thing to integrate"
    
    position = y[0:3]
    velocity = y[3:6]
    
    dydt = np.empty(6)
    
    dydt[0:3] = velocity
    dydt[3:6] = np.array([0, 0, -9.81]) #force goes here...
    # Don't forget to divide by m
    
    return dydt 

m = 1
theta = np.pi/3
t = np.linspace(0, 1)
y0 = np.array([0, 0, 0, 5*np.cos(theta), 0, 5*np.sin(theta)])
```


```python
y = odeint(main, y0, t)
plt.plot(t, y[:, 2])
```


```python
def df(f, x):
    h = 0.0001
    return (f(x + h) - f(x)/h) 

# Whenever possible, use polyderiv if possible
```


```python
# If number of counts is n, then the uncertainty is sqrt(n) cough Poisson distribution
# Error of the mean is always stdev/sqrt(n) by CLT
```

## Curve Fitting Syntax

To use curve fit please pass in a function with params


```python
import numpy as np 
from scipy.optimize import curve_fit

x = np.array([1, 2, 3, 9])
y = np.array([1, 4, 1, 3])

def fit_func(x, a, b):
    return a*x + b

popt = curve_fit(fit_func, x, y)
popt[0]

plt.scatter(x, y)

x_cont = np.linspace(1, 9)
plt.plot(x_cont, fit_func(x_cont, *popt[0]))
```

## Minimization Syntax

To use minimize must pass in a loss function. Note we already plug in the data


```python
import numpy as np 
from scipy.optimize import minimize

x = np.array([1, 2, 3, 9])
y = np.array([1, 4, 1, 3])

def loss_func(params, x, y):  
    a, b = params
    
    fit_func = lambda x: a*x + b
    return np.sum((fit_func(x) - y)**2)
```


```python
params0 = np.array([1, 1])
popt_minimization = minimize(loss_func, params0, args=(x, y))
```


```python
plt.scatter(x, y)
plt.plot(x_cont, fit_func(x_cont, *popt_minimization.x))
```

## Least Squares Syntax


```python
import numpy as np 
from scipy.optimize import least_squares

def res_func(params, x, y):
    a, b = params
    fit_func = lambda x: a*x + b

    return (fit_func(x) - y)

params0 = np.array([1, 1])
popt_minimization = least_squares(res_func, params0, args=(x, y))
popt_minimization

plt.scatter(x, y)
plt.plot(x_cont, fit_func(x_cont, *popt_minimization.x))
```

## Polyfit Syntax


```python
import numpy as np 

x = np.array([1, 2, 3, 9])
y = np.array([1, 4, 1, 3])

plt.scatter(x, y)

x_cont = np.linspace(1, 9)

popt = np.polyfit(x, y, 1)
popt

plt.plot(x_cont, np.poly1d(popt)(x_cont))

# Poly fit is the easiest to use
```

## Testing


```python

    
```


```python
myfunc()
```

    3



```python
a = np.arange(4, 10)
wh = np.where(a > 7)
a[wh]
```




    array([8, 9])




```python
a[a > 7]
```




    array([8, 9])




```python
b = np.arange(10, 30).reshape(-1, 4)
b
```




    array([[10, 11, 12, 13],
           [14, 15, 16, 17],
           [18, 19, 20, 21],
           [22, 23, 24, 25],
           [26, 27, 28, 29]])




```python
wh = np.where(b > 15)
b[wh[0], wh[1]]
```




    array([16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29])




```python
(b > 15)

```




    array([[False, False, False, False],
           [False, False,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True]])




```python
(b < 17) & (b > 15)
# & is a ufunc
# and comparing the whole array

```




    array([[False, False, False, False],
           [False, False,  True, False],
           [False, False, False, False],
           [False, False, False, False],
           [False, False, False, False]])




```python
xdata = np.linspace(1, 10, 10)
ydata = np.exp(xdata)

x_cont = np.linspace(1, 10)

plt.plot(x_cont, interp1d(xdata, ydata, kind='linear')(x_cont))
plt.plot(x_cont, interp1d(xdata, ydata, kind='cubic')(x_cont))
# sometimes need to turn on extrapolate if minimize is complaing
```


```python

```

d1fa1acd3c693b2b8d7aded2379cf7c4f8f6f853
ipynb
Jupyter Notebook
Assignments/Useful Functions.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Useful Functions.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Assignments/Useful Functions.ipynb
hanzhihua72/phys-420
748d29b55d57680212b15bb70879a24b79cb16a9
MIT
Qwen/Qwen-72B
1. YES
2. YES