Read Codify Learn

 # Q 7A import numpy as np from scipy.special import legendre as P theta = np.pi n= 2 lhs = np.sin((n+1)*theta)/np.sin(theta) rhs = sum([np.poly1d(P(l))(np.cos(theta)) * np.poly1d(P(n-l))(np.cos(theta)) for l in range(0, n+1)]) print("LHS = ",lhs) print("RHS = ",rhs) print("Difference = ", abs(lhs-rhs))

 # Q 6A # Integral e^-x2 from 0 --> infinity import numpy as np from scipy.integrate import quad, simps import matplotlib.pyplot as plt f = lambda x: np.exp(-x**2) res = quad(f, 0, np.inf) print("Value of integral is = ", round(res[0], 4))

 # Q 5B import numpy as np from scipy.integrate import quad, simps from scipy.signal import square, sawtooth import matplotlib.pyplot as plt x = np.linspace(-np.pi, np.pi, 1001) f1 = abs(np.sin(x)) L = 2*np.pi # Evaluation of Fourier coefficients a0 = (2.0/L)*simps(f1, x) an = lambda n: (2.0/L) * simps(f1*np.cos(n*x), x) bn = lambda n: (2.0/L) * simps(f1*np.sin(n*x), x) # Summing the serier R = 10 xn = np.linspace(-R, R, 1001) s = 0.5*a0 +  sum([an(n)*np.cos(n*xn) + bn(n)*np.sin(n*xn) for n in range(1, 50)]) signal = abs(np.sin(xn)) #Plotting the Functions plt.plot(xn, signal,'r-' ,label="Original Signal") plt.plot(xn, s,'ko',ms=2, label = "Fourier Approximation") plt.legend() plt.show()

 import numpy as np from scipy.integrate import quad import matplotlib.pyplot as plt f = lambda x: np.sin(x**2) A = [] x = np.arange(0.01, 10, 0.01) for i in x:          res = quad(f, 0, i)     A.append(res[0]) print(sum(A)) plt.plot(x[:600], A[:600],label="Fresnel Integral Si(x)") plt.grid() plt.legend() plt.show() print("From the plot, we find that the solution slowly converges to 0.5")

 # Q 4B # Question is nothing but 1-D Heat Equation with heating at middle # Solve by using finite difference method import numpy as np import matplotlib.pyplot as plt from scipy.integrate import quad x = np.linspace(0, 1, 101) u = np.zeros(101) u[0] = u[100] = 0       # Boundary Conditions u[50] = 1               # Initial Condition for t in range(100):     for i in range(1,100):         u[i] += (u[i+1] + u[i-1] -2*u[i])/4 plt.plot(x, u, 'r-', lw=3) plt.xlabel("Length of rod") plt.ylabel("Temperature (in reduced units)") plt.title("Temperature Profile of Rod") plt.show()

 # Q 4A # Fourier Series Sawtooth import numpy as np from scipy.signal import sawtooth from scipy.integrate import simps import matplotlib.pyplot as plt t = np.linspace(0,1, 1001) f1 = sawtooth(2*np.pi*2*t,1) T = 0.5 #Calculation Of Coefficients a0 = (2./T) * simps(f1, t) an = lambda n : (2./T) * simps(f1*np.cos(2*np.pi*2*n*t), t) bn = lambda n : (2./T) * simps(f1*np.sin(2*np.pi*2*n*t), t) tn = np.linspace(0, 6, 1001) s =0.5*a0 + sum([an(n)*np.cos(2*np.pi*2*n*tn) + bn(n)*np.sin(2*np.pi*2*n*tn) for n in range(10)]) plt.plot(tn, s, 'r',label="Fourier approximation") plt.plot(tn, 2*sawtooth(2*np.pi*2*tn,1) ,label= "Original Signal") plt.legend() plt.show()

 # Q 3B import numpy as np from scipy.special import legendre from scipy.integrate import quad # 1. for m = n (= 6,say) m=6 n=6 f1 = legendre(m) * legendre(n) y = quad(f1, -1, 1) print("\nIntegral Pm * Pn for m = n = 6 is:",y[0])        # Verifying that integral of P(m)*P(n) from -infinity to +infinity for m=n equals 2/(2m+1) print("2/(2m + 1); m = 6 is = " ,2/(2*m + 1)) # 2. for  m  !=  n m = 3 n = 2 f2 = legendre(m) * legendre(n) y2 = quad(f2, -1, 1)                                # Verifying that integral of P(m)*P(n) from -infinity to +infinity for m!=n equals 0 print("\nIntegral Pm * Pn for m not = n is:",y2[0])

 # Q 3A # d2x/dt2 + gamma * dx/dt + kx = Fcos wt # Exact Sol of Amplitude = F / sqrt((k-w)^2 + gamma^2*w^2) # take gamma = 0.5 , k = 5 , F = 1 import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint k  = 5 gam = 0.5 F = 1. def forced(x ,t):     return [x[1], F*np.cos(w*t) - gam*x[1] - k*x[0]] x0 = [0., 1.] t = np.linspace(0, 100, 1001) freq = np.linspace(0.5, 5, 1001) A = [] for w in freq:     sol = odeint(forced, x0, t)     X = sol[:,0]     amp = np.max(X[-200:])     A.append(amp) exact = F/(np.sqrt((k-freq**2)**2 + gam**2*freq**2))      plt.plot(freq[::2], A[::2], '--', label="soultion by odeint") plt.plot(freq, exact, 'o', ms=0.5, label ="exact solution") plt.xlabel("frequency") plt.ylabel("Amplitude") plt.legend() plt.show()

 # Q 2B import numpy as np from math import factorial from scipy.special import hermite as H from scipy.integrate import quad import matplotlib.pyplot as plt x =np.linspace(-5,5, 1001) def psi(n):     return 1/np.sqrt(2**n * factorial(n)) * (1/np.pi)**(1/4) * np.exp(-x**2/2) * H(n)(x) fig, ax = plt.subplots(1,3, sharex=True) ax[0].plot(x, psi(2),'r') ax[0].set_title("$\psi_2$") ax[0].axhline() ax[1].plot(x, psi(1),'b') ax[1].set_title("$\psi_1$") ax[1].axhline() ax[2].plot(x, psi(0),'orange') ax[2].set_title("$\psi_0$") ax[2].axhline() plt.show()

 # Q 2A #d2x/dt2 + gamma*dx/dt + k*x = 0 #take gamma = 0.4, k = 1 and solve using suitable scipy function import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt gam = 0.4 k =  1. def damp(x, t):     return [x[1], -k*x[0]-gam*x[1]] x0 = [1.0, 0.0] t = np.linspace(0,20,1001) sol = odeint(damp, x0, t) plt.plot(t, sol[:,0], 'r', label="x") plt.xlabel("time(t)") plt.ylabel("x") plt.axhline() plt.title("Damped Harmonic Oscillator") plt.legend() plt.show()

 # Q 1B import numpy as np import matplotlib.pyplot as plt from scipy.integrate import simps a = 1 b = 2 F = lambda t: np.exp(-a*t**2) F = np.vectorize(F) G = lambda t: np.exp(-b*t**2) G = np.vectorize(G) conv = [] t = np.linspace(-5, 5, 1000) for i in t:     t1 = np.linspace(-5.001, i, 1000)     h = F(t1) * G(i - t1)     conv.append(simps(h,t1)) plt.plot(t, F(t), label = "F(t)") plt.plot(t, G(t), label = "G(t)") plt.plot(t, conv, label = "Convolution") plt.title("Ques 1B : Verify that convolution of 2 Gaussians is also a Gaussian") plt.legend() plt.savefig("1B.png", dpi = 600) plt.show()

#Q 1A import numpy as np from scipy.integrate import quad import matplotlib.pyplot as plt # Defining the Gaussian Function def Gaussian(mu, sigma, x):     return (1./(np.sqrt(2*np.pi)*sigma)) * np.exp(-(x-mu)**2/(2*sigma**2)) x = np.linspace(-5,5,1001) mu = 0 sigma = 0.5 # Plotting the Gaussian Function plt.plot(x, Gaussian(mu,sigma,x), label=''r'$\frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$') plt.xlabel('x') plt.ylabel('f(x)') plt.title("Ques. 1A Gaussian Function") plt.legend() plt.savefig("1A.png", dpi=600) plt.show() # Evaluating the Integral mu = 1 sigma = 0.5 def f(x):     return (1./(np.sqrt(2*np.pi)*sigma)) * np.exp(-(x-mu)**2/(2*sigma**2)) x = np.linspace(-5,5,1001) result_1 = quad(f, -np.inf, np.inf) result_2 = 1./(sigma * np.sqrt(2*np.pi)) * result_1[0] print("Value of the integral is = ", result_2)