<jupyter_start><jupyter_text>SciPy Tutorial<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import scipy as sp<jupyter_output><empty_output><jupyter_text>Basic  Optimization<jupyter_code>from scipy.optimize import minimize<jupyter_output><empty_output><jupyter_text>Minimize $f(x) = (x-3)^2$<jupyter_code>def f(x):
    return (x-3)**2    # defining the function

res = minimize(f, 1)    # giving an initial guess, can be any value
res.x    # gives an array
res.x[0]   # taking the value of x<jupyter_output><empty_output><jupyter_text>Minimize $f(x, y) = (x-1)^2 + (y-2.5)^2$ subject to:$$x-2y+2 \geq 0 $$$$-x-2y+6 \geq 0$$$$-x + 2y + 2 \geq 0$$$$x \geq 0 $$$$y \geq 0 $$![opt_reg](img/optimization/opt_reg.png)* 2D function takes in vector $x$* Constraints must be specified as $g_i(x) \geq 0$* Bounds specified as rectangular<jupyter_code># defining the equation in a variable
f = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2    # x[0], x[1] represent x and y axis respectively

# giving the constraints in a tuple of dictionary
cons =  ({'type': 'ineq', 'fun': lambda x: x[0] - 2*x[1] + 2},
                {'type': 'ineq', 'fun': lambda x: -x[0] - 2*x[1] + 6},
                {'type': 'ineq', 'fun': lambda x: -x[0] + 2*x[1] + 2})

# giving the boundaries where None means infinity
bnds = ((0, None), (0, None))

res = minimize (f, (1, 0), bounds = bnds, constraints = cons)

res.x
# diff methods to minimize
minimize?<jupyter_output>[1;31mSignature:[0m
[0mminimize[0m[1;33m([0m[1;33m
[0m    [0mfun[0m[1;33m,[0m[1;33m
[0m    [0mx0[0m[1;33m,[0m[1;33m
[0m    [0margs[0m[1;33m=[0m[1;33m([0m[1;33m)[0m[1;33m,[0m[1;33m
[0m    [0mmethod[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mjac[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mhess[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mhessp[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mbounds[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mconstraints[0m[1;33m=[0m[1;33m([0m[1;33m)[0m[1;33m,[0m[1;33m
[0m    [0mtol[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0mcallback[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m    [0moptions[0m[1;33m=[0m[1;32mNone[0m[1;33m,[0m[1;33m
[0m[1;33m)[0m[1;33m[0m[1;33m[0m[0m
[1;31mDocstring:[0m
Minimization of scalar function of one or more variables.

Parameters
----------
fun : callable[...]<jupyter_text>![minimize_methods](img/optimization/minimize_methods.svg)<jupyter_code>x = np.linspace(0, 10, 10)
y = x**2 * np.sin(x)
plt.scatter(x,y)<jupyter_output><empty_output><jupyter_text>Now, I want to know the values in between.<jupyter_code>from scipy.interpolate import interp1d    # importing the module to interpolate (joining the points)
f = interp1d(x, y, kind='cubic')    # we can give linear also resulting joining the points using straight lines
x_dense = np.linspace(0, 10, 100)   # there will be total of 100 points in between the 10 points. more points, better interpolation
y_dense = f(x_dense)   # will take corresponding value of x
plt.plot(x_dense, y_dense)
plt.scatter(x, y)<jupyter_output><empty_output><jupyter_text>Curve Fitting*Trivial Example*<jupyter_code>x_data = np.linspace(0, 10, 10)
y_data = 3*x_data**2 + 2
plt.scatter(x_data, y_data)<jupyter_output><empty_output><jupyter_text>Want to fit the data to the curve $y = ax^2 + b$. The main goal here is determine the values of $a$ and $b$<jupyter_code># importing the method
from scipy.optimize import curve_fit

# defining the function
def func(x, a, b):
    return a*x**2 + b


popt, pcov = curve_fit(func, x_data, y_data, p0=(1, 1))

popt<jupyter_output><empty_output><jupyter_text>*Real Example from Physics* The equation for spring motion is $y(t) = A\cos(\omega t+\phi)$. Want to find the natural frequency of oscillation $\omega$ for the spring. You collect the data<jupyter_code>t_data = np.array([0.,  0.34482759,  0.68965517,  1.03448276,  1.37931034,
                   1.72413793,  2.06896552,  2.4137931,  2.75862069,  3.10344828,
                   3.44827586,  3.79310345,  4.13793103,  4.48275862,  4.82758621,
                   5.17241379,  5.51724138,  5.86206897,  6.20689655,  6.55172414,
                   6.89655172,  7.24137931,  7.5862069,  7.93103448,  8.27586207,
                   8.62068966,  8.96551724,  9.31034483,  9.65517241, 10.])
y_data = np.array([4.3303953,  1.61137995, -2.15418696, -3.90137249, -1.67259042,
                   2.16884383,  3.86635998,  1.85194506, -1.8489224, -3.96560495,
                   -2.13385255,  1.59425817,  4.06145238,  1.89300594, -1.76870297,
                   -4.26791226, -2.46874133,  1.37019912,  4.24945607,  2.27038039,
                   -1.50299303, -3.46774049, -2.50845488,  1.20022052,  3.81633703,
                   2.91511556, -1.24569189, -3.72716214, -2.54549857,  0.87262548])

plt.plot(t_data, y_data, 'o--')<jupyter_output><empty_output><jupyter_text>$\omega = 2 \pi f$, $f = 1/T$, and $T \approx 2$ seconds. Thus good initial guess is*   $\omega = 2 \pi(1/2) = \pi$*   $A = 4$*   $\phi = 0$<jupyter_code>from scipy.optimize import curve_fit


def func(x, A, w, phi):
    return A*np.cos(w*x+phi)


# Guess Should Be Logical, Otherwise The Fitting Will Be Horrible
popt, pcov = curve_fit(func, t_data, y_data, p0=(4, np.pi, 0))

popt   # gives value most probable value of A, w and phi respectively

#  parameters
A, w, phi = popt
t = np.linspace(0, 10, 100)
y = func(t, A, w, phi)
plt.scatter(t_data, y_data)
plt.plot(t, y)
np.diag(pcov)   # the deviation of parameters
np.sqrt(np.diag(pcov))   # error in calculation of parameters<jupyter_output><empty_output><jupyter_text>Special Functions **Legendre Polynomials:**  $P_l(x)$: satisfy $(1-x^2)y'' -2xy'+l(l+1)y = 0$*   The angular component of the spherical Schrodinger Equation which permits non-infinite<jupyter_code># importing and defining the legendre polynomial equation
from scipy.special import legendre
x = np.linspace(0, 1, 100)
plt.plot(x, legendre(6)(x))
plt.show()
legendre(6)<jupyter_output><empty_output><jupyter_text>**Bessel functions:** $J_{\alpha}(x)$: satisfy $x^2 y'' + xy' + (x^2-\alpha^2)y = 0$*   Laplace's Equation in Polar Coordinates<jupyter_code>from scipy.special import jv
x = np.linspace(0, 10, 100)
plt.plot(x, jv(3,x))
plt.show()<jupyter_output><empty_output><jupyter_text>Calculus  Differentiation Differentiation of the following equation -  $$x^2\sin{(2x)}e^{-x}$$<jupyter_code>from scipy.misc import derivative
def f(x):
    return x**2 * np.sin(2*x) * np.exp(-x)

x = np.linspace(0, 1, 100)
plt.plot(x, f(x))   # plotting the function
plt.plot(x, derivative(f, x, dx=1e-6))   # 1e-6 separation between two points
plt.plot(x, derivative(f, x, dx=1e-6, n=2))   # 2nd order derivative<jupyter_output><empty_output><jupyter_text>Integration Integration of the following single integral- $$\int_{0}^1 x^2\sin{(2x)}e^{-x}dx$$<jupyter_code>from scipy.integrate import quad   # quad method handels the single integration part
def integrand(x): return x**2 * np.sin(2*x) * np.exp(-x)
integral, integral_error = quad(integrand, 0, 1)

# Following lines mean same as previous (using lambda function)

# integrand = lambda x: x**2 * np.sin(2*x) * np.exp(-x)
# integral, integral_error = quad(integrand, 0, 1)

integral   # value of the integration
integral_error   # error in the integration value calculation<jupyter_output><empty_output><jupyter_text>Integration of the following double integral- $$\int_{0}^1\int_{-x}^{x^2}\sin{(x+y^2)}\,dydx$$<jupyter_code>from scipy.integrate import dblquad   # dblquad method handels the double integration part
def integrand(y, x): return np.sin(x+y**2)
def lwr_y(x): return -x
def upr_y(x): return x**2
integral, integral_error = dblquad(integrand, 0, 1, lwr_y, upr_y)

# Following lines mean same as previous (using lambda function)

# integrand = lambda y, x: np.sin(x+y**2)
# lwr_y = lambda x: -x
# upr_y = lambda x: x**2
# integral, integral_error = dblquad(integrand, 0, 1, lwr_y, upr_y)
integral
integral_error<jupyter_output><empty_output><jupyter_text>Differential Equations  Single First order ODEs**ODE** $\equiv$ **Ordinary Differential Equations**Air friction while falling $$v^{'}-\alpha v^2+ \beta = 0 \hspace {10mm} v(0) = 0$$<jupyter_code>from scipy.integrate import odeint   # importing the library for solving ordinary differential equation
def dvdt(v, t):
    return 3*v**2 - 5 

v0 = 0   # initial condition<jupyter_output><empty_output><jupyter_text>I'm returning the value of $\alpha$ and $\beta$, 3 and 5 respectively (arbitary) in the equation $v^{'}=\alpha v^2+\beta$<jupyter_code>t = np.linspace(0, 1, 100)   # within a timerange (unitless in SciPy, only value we get) of 0 and 1
sol = odeint(dvdt, v0, t)
sol.T[0]   # values for the solution after taking the transpose (changing into a single row matrix)
plt.plot(t, sol.T[0])   # plot for the solution<jupyter_output><empty_output><jupyter_text>Coupled First Order ODEs $$y^{'}=y_1+y^2_2+3x\hspace{10mm}y_1(0)=0$$$$y^{'}_2=3y_1+y^3_2-\cos(x)\hspace{10mm}y_2(0)=0$$Letting $S=(y_1, y_2)$ we need to write a function that returns $dS/dx = (dy_1/dx, dy_2/dx)$. The function $dS/dx$ can take in $S=(y_1, y_2)$ and $x$<jupyter_code>def dSdx(S, x):
    y1, y2 = S
    return [y1 + y2**2 + 3*x, 3*y1 + y2**3 - np.cos(x)]   # returns (dy1/dx and dy2/dx)

# initial conditions
y1_0 = 0
y2_0 = 0

# giving S the value of initial conditions
S_0 = (y1_0, y2_0)
x = np.linspace(0, 1, 100)   # solvind for 100 values between the value of x 0 and 1
sol = odeint(dSdx, S_0, x)
sol   # it gives an array of values of y1 and y2 at every value of x between 0 and 1
y1 = sol.T[0]
y2 = sol.T[1]
# gives the plot for above coupled first order ODE
plt.plot(x, y1)
plt.plot(x, y2)<jupyter_output><empty_output><jupyter_text>Second Order ODEEquation for a pendulum$$\theta'' - \sin(\theta) = 0$$* Scipy can only solve coupled first order ODEs, but any second order ODE can be turned into two coupled first order ODEs. The same thing goes for higher order ODEs.Define $\omega = d\theta/dt$ so that one has the following coupled ODEs$$d \omega / dt = \sin(\theta)$$ $$d \theta / dt = \omega $$Let $S = (\theta, \omega)$<jupyter_code>def dSdt(S, t):
    theta, omega = S
    return [omega, np.sin(theta)]

theta0 = np.pi/4
omega0 = 0
S0 = (theta0, omega0)

t = np.linspace(0, 20, 100)
sol = odeint(dSdt, S0, t)
theta, omega = sol.T

plt.plot(t, theta)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16669841874642788, 'fraction_numerical': 0.09820918270146695, 'mean_word_length': 4.491108786610878, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 76, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '597212', 'license': 'MIT', 'n_tokens_mistral': 5006, 'n_tokens_neox': 4341, 'n_words': 1041, 'path': 'SciPy/SciPy Tutorial.ipynb', 'repo_name': 'aritraroy24/Learning_NumPy_SciPy', 'size': '10498'}