API tools faq
paste
Advertisement
yusufbrima

Python Assignment

yusufbrima
Oct 18th, 2020
5,285
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 11.64 KB | None | 0 0
raw download clone embed print report
  1. #!/usr/bin/env python
  2. # coding: utf-8
  3.  
  4. # In[31]:
  5.  
  6.  
  7. get_ipython().run_cell_magic('javascript', '', 'MathJax.Hub.Config({\n    TeX: { equationNumbers: { autoNumber: "AMS" } }\n});')
  8.  
  9.  
  10. # # $$\text{Exercise 1}$$
  11. # 1. Let us consider the array $M$ defined below
  12.  
  13. # In[47]:
  14.  
  15.  
  16. import numpy as np
  17. import matplotlib.pyplot as plt
  18. M = np.array([[ 80842, 333008, 202553, 140037,  81969,  63857,  42105, 261540],
  19.        [481981, 176739, 489984, 326386, 110795, 394863,  25024,  38317],
  20.        [ 49982, 408830, 485118,  16119, 407675, 231729, 265455, 109413],
  21.        [103399, 174677, 343356, 301717, 224120, 401101, 140473, 254634],
  22.        [112262,  25063, 108262, 375059, 406983, 208947, 115641, 296685],
  23.        [444899, 129585, 171318, 313094, 425041, 188411, 335140, 141681],
  24.        [ 59641, 211420, 287650,   8973, 477425, 382803, 465168,   3975],
  25.        [ 32213, 160603, 275485, 388234, 246225,  56174, 244097,   9350],
  26.        [496966, 225516, 273338,  73335, 283013, 212813,  38175, 282399],
  27.        [318413, 337639, 379802, 198049, 101115, 419547, 260219, 325793],
  28.        [148593, 425024, 348570, 117968, 107007,  52547, 180346, 178760],
  29.        [305186, 262153,  11835, 449971, 494184, 472031, 353049, 476442],
  30.        [ 35455, 191553, 384154,  29917, 187599,  68912, 428968,  69638],
  31.        [ 20140, 220691, 163112,  47862, 474879,  15713, 279173, 161443],
  32.        [ 76851, 450401, 193338, 189948, 211826, 384657,  17636, 189954],
  33.        [182580, 403458,  63004,  87202,  17143, 308769, 316192, 408259],
  34.        [379223, 322891, 302814, 309594,  13650,  68095, 327748,  27138],
  35.        [ 33346, 202960, 221234,  61986,  96524, 230606, 266487, 299766],
  36.        [ 46264,  33510, 434647,  26750, 412754, 346058, 263454, 270444],
  37.        [182770, 486269, 402948, 105055, 463911, 378799, 329053,  16793],
  38.        [379586, 293162, 121602, 144605, 105553,  71148,  40126, 183534],
  39.        [493338,   5241, 498084, 189451, 219466,  12869, 482845,  16529],
  40.        [145955, 461866, 427009,  70841, 318111, 199631,  99330, 407887],
  41.        [133826, 214053, 295543, 339707, 199545, 194494, 304822, 337793],
  42.        [395524, 206326, 308488, 297633, 472409, 487606, 422244, 443487],
  43.        [464749,  75275, 299486, 248984,  76076, 197544,  75103, 394086],
  44.        [404365,  67309, 267556,  96483, 166099, 203867,  57120, 211249],
  45.        [ 86992, 286426, 130451, 419955, 109839, 332663, 101095,  50820],
  46.        [108699, 436141, 333670,  45890, 425798,   6967, 155851, 360687],
  47.        [252907, 121138,  81509, 367876, 471831, 226950, 349068, 197092]])
  48. M
  49.  
  50.  
  51. # 1. What is the average value of the second and fourth columns?
  52.  
  53. # In[38]:
  54.  
  55.  
  56. np.average(M[:,1:4:2])
  57.  
  58.  
  59. # In[8]:
  60.  
  61.  
  62. np.average(M[:,1:4:2])
  63.  
  64.  
  65. # 2. What is the average value of the last 5 rows of the second and fourth columns?
  66.  
  67. # In[10]:
  68.  
  69.  
  70. np.average(M[-5:,1:4:2])
  71.  
  72.  
  73. # 3. Update $M$ by Replacing all its even integers by 0.
  74.  
  75. # In[11]:
  76.  
  77.  
  78. M[M%2==0]=0
  79. M
  80.  
  81.  
  82. # # $$\text{Exercise 2}$$
  83. #
  84. # Let $\{ x_k\}$ be a partition of $[a,b]$ such that $a=x_0<x_1<\cdots<x_{N-1}<x_N=b$ and $H$ be the length of the $k$-th subinterval ($H = x_k - x_{k-1}$),
  85. # then we have
  86. # $$\int_a^bf(x)dx \approx \sum_{k=1}^N \frac{f(x_{k-1})+f(x_k)}{2}H = A$$
  87. #
  88. #
  89. # Write a function named <b>Trap</b> that takes $a,b,N, f$ as inputs and return A
  90.  
  91. # In[39]:
  92.  
  93.  
  94. import numpy as np
  95. def Trap(a,b,N,f):
  96.     x = np.linspace(a,b,N+1)
  97.     H = x[1]-x[0]
  98.     A = 0
  99.     for k in range (1,N+1):
  100.         A+=(H/2)*(f(x[k-1])+f(x[k]))
  101.     return A
  102. Trap(0,10,20,lambda x:x**2)
  103.  
  104.  
  105. # # $$\text{Exercise 3}$$
  106. #
  107. # Let us consider the function $f$ defined on the set of integer as follows
  108. #
  109. # \begin{equation}
  110. # f(n) =
  111. #   \begin{cases}
  112. #     n/2       & \quad \text{if } n \text{ is even}\\
  113. #     -(n)/3 +1  & \quad \text{if } n \text{ is odd and divisible by 3}\\
  114. #     (n+1)/2 +1 & \quad \text{else}
  115. #   \end{cases}
  116. # \end{equation}
  117. #
  118. # Write a python function named <b> ComptFunc </b> that takes two intergers $a$ and $b$ , $a<b$ and return  an array $flist$ of all  $f(n), n \in nlist$ where $nlist$ is the list of all the integers in the interval $[a,b]$.
  119.  
  120. # In[40]:
  121.  
  122.  
  123. def f(n):
  124.     if n%2==0:
  125.         return n/2
  126.     elif n%2==1 and n%3==0:
  127.         return (-(n)/3)+1
  128.     else:
  129.         return ((n+1)/2)+1
  130. vf = np.vectorize(f)
  131. def CompFunc(a,b):
  132.     if a<b:
  133.         nlist = np.arange(a,b+1)
  134.         flist = vf(nlist)
  135.         return flist
  136. CompFunc(0,9)
  137.  
  138.  
  139. # # $$\text{Exercise 4}$$
  140. #
  141. # Write a python code to create the following numpy arrays
  142. # $$
  143. # A = \begin{pmatrix}
  144. # 1 & 4 & 6 \\
  145. # 0 & -3 & 2 \\
  146. # -2 & -2 & -2 
  147. # \end{pmatrix}, \quad
  148. # B = \begin{pmatrix}
  149. # 2 & -1 & 0 \\
  150. # 2 & -1 & 0 \\
  151. # 2 & -3 & 1
  152. # \end{pmatrix},
  153. # $$
  154. # and compute
  155. # - $A-B$,
  156. # - $4A + B$,
  157. # - $trace(A)$,
  158. # - $B^t$ the transpose of $B$,
  159. # - $AB$,
  160. # - $BA^t$,
  161. # - the determinant of $A$
  162.  
  163. # In[13]:
  164.  
  165.  
  166. A = np.array([[1,4,6],[0,-3,2],[-2,-2,-2]])
  167. B = np.array([[2,-1,0],[2,-1,0],[2,-3,1]])
  168.  
  169.  
  170. # In[14]:
  171.  
  172.  
  173. A-B
  174.  
  175.  
  176. # In[15]:
  177.  
  178.  
  179. 4*A+B
  180.  
  181.  
  182. # In[16]:
  183.  
  184.  
  185. np.trace(A)
  186.  
  187.  
  188. # In[17]:
  189.  
  190.  
  191. B.T
  192.  
  193.  
  194. # In[ ]:
  195.  
  196.  
  197. A
  198.  
  199.  
  200. # In[ ]:
  201.  
  202.  
  203.  
  204.  
  205.  
  206. # # $$\text{Exercise 5}$$
  207. #
  208. # write a python code to solve the system of equation  
  209. # \begin{equation}\label{sysEqu}
  210. # \left\lbrace
  211. # \begin{array}{ll}
  212. # 2x + y -2z &= 3,\\
  213. # x-y-z &= 0,\\
  214. # x+y+3z &=12
  215. # \end{array}\right.
  216. # \end{equation}
  217.  
  218. # In[42]:
  219.  
  220.  
  221. A = np.array([[2,1,-2],[1,-1,-1],[1,1,3]])
  222. B = np.array([3,0,12])
  223. P = np.linalg.solve(A,B)
  224. P
  225.  
  226.  
  227. # # $$\text{Exercise 6}$$
  228. #
  229. # Let us consider the sequence $U_n$ given by
  230. # \begin{equation}\label{fib}
  231. # \left\lbrace
  232. # \begin{array}{ll}
  233. # U_0 &= 1,\\
  234. # U_1 &= 2,\\
  235. # U_{n} &=-3U_{n-1} +U_{n-2}, \;\; \forall\; n=2,3,4\cdots
  236. # \end{array}\right.
  237. # \end{equation}
  238. #
  239. # Write a python function named  <b>SeqTerms</b> that takes as input an integer $n,\;\;n\geq 0$ and return an array of the first $n+1$ terms (i.e. $U_0, \cdots, U_{n}$) of the sequence \eqref{fib}.
  240.  
  241. # In[118]:
  242.  
  243.  
  244. def SeqTerms(n):
  245.     U =  np.array([])
  246.     if(n >= 0):
  247.             if(n == 0):
  248.                 return U
  249.             elif(n == 1):
  250.                 return np.append(U,1)
  251.             elif(n == 2):
  252.                 return np.append(U,[1,2])
  253.             else:
  254.                 U = np.append(U,[1,2])
  255.                 for i in range(2,n):
  256.                     U = np.append(U,-3*U[i-1] + U[i-2])
  257.                 return U.astype(int)
  258. SeqTerms(10)
  259.  
  260.  
  261.  
  262.  
  263. # In[71]:
  264.  
  265.  
  266. len([])
  267.  
  268.  
  269. # # $$\text{Exercise 7}$$
  270. #
  271. # Let $\{ x_k\}$ be a partition of $[0,1]$ such that $0=x_0<x_1<\cdots<x_{99}<x_{100}=1$ and $H$ be the length of the $k$-th subinterval ($H = x_k - x_{k-1}$).
  272. #
  273. # 1. Write a python code that use Euler’s method to solve the initial value problem
  274. # \begin{equation}\label{eul1}
  275. #   \begin{cases}
  276. #     y' = \dfrac{2-2xy}{x^2+1},      & \quad \text{on } [0, 1]\\\\
  277. #     y(0) = 1,
  278. #   \end{cases}
  279. # \end{equation}
  280. # i.e. generates array of  $y_k\approx g(x_k)=g_k$ where $g$ is the exact solution of \eqref{eul1}
  281. #
  282. # 2. The exact solution of the initial value problem is given by
  283. # \begin{equation}\label{exact}
  284. # g(x) = \dfrac{2x+1}{x^2+1}
  285. # \end{equation}
  286. #
  287. # use $subplot$ to plot side by side
  288. # - the exact and appromixate solution in the same window (i.e. ($y_k$  vs $x_k$ and $g_k$ vs $x_k$))
  289. # - the error $e_k = |y_k -g_k| $ against $x_k$.
  290.  
  291. # In[108]:
  292.  
  293.  
  294. f = lambda x,y: (2-(2*x*y))/(x**2 +1)
  295. g = lambda x: (2*x +1)/(x**2 +1)
  296. z = lambda f,g: abs(f-g)
  297. def Euler(a,b,n,f):
  298.     x = np.linspace(a,b,n)
  299.     H = (b-a)/n
  300.     y = np.zeros(n)
  301.     y[0] = 1
  302.     for i in range(1,n):
  303.         y[i] = y[i-1] + (H * f(x[i-1],y[i-1]))
  304.     return x,y
  305.  
  306.  
  307. # In[135]:
  308.  
  309.  
  310. x,y = Euler(0,1,100,f) #approximate solution
  311. plt.plot(x,y)
  312. y_0 = g(x) #actual solution
  313.  
  314. #plt.plot(x,y_0)
  315. #plt.legend(["f(x)","g(x)"])
  316. #plt.title("Euler's Method of function approximation")
  317. #plt.xlabel("x")
  318. #plt.ylabel("f(x),g(x)")
  319.  
  320. err = z(y,y_0)
  321. #print(err)
  322.  
  323.  
  324. fig, (ax1, ax2) = plt.subplots(1, 2)
  325. fig.suptitle('Function approximation')
  326. ax1.plot(x, y)
  327. ax1.title.set_text('Euler\'s Method of function approximation')
  328. ax2.plot(x, err)
  329. ax2.title.set_text('Approximation Error')
  330.  
  331.  
  332. # In[70]:
  333.  
  334.  
  335. x = np.linspace(0,10,100)
  336.  
  337. y = np.exp(x)
  338. z = np.exp(-x) + np.sin(x)
  339.  
  340. fig, (ax1, ax2) = plt.subplots(1, 2)
  341. fig.suptitle('Horizontally stacked subplots')
  342. ax1.plot(x, y)
  343. ax1.title.set_text('First Plot')
  344. ax2.plot(x, z)
  345. ax2.title.set_text('Second Plot')
  346.  
  347.  
  348. # # $$\text{Exercise 8}\quad (\text{generalization of previous exercise})$$
  349. #
  350. # Let $\{ x_k\}$ be a partition of $[a,b]$ such that $a=x_0<x_1<\cdots<x_{99}<x_{N}=b$ and $H$ be the constant length of the $k$-th subinterval ($H = x_k - x_{k-1}$). Let us consider initial value problem
  351. #
  352. # \begin{equation}\label{eul2}
  353. #   \begin{cases}
  354. #     y' = f(x,y),      & \quad \text{on } [a, b]\\\\
  355. #     y(a) = c,
  356. #   \end{cases}
  357. # \end{equation}
  358. #
  359. # 1. Write a python function <b> EulerMethod </b> that takes $a,b,c,N,$ and $f$ and plot the approximate solution of \eqref{eul2} using Euler method.
  360.  
  361. # In[116]:
  362.  
  363.  
  364. def EulerMethod(a,b,c,N,f):
  365.     x = np.linspace(a,b,N)
  366.     H = (b-a)/N
  367.     y = np.zeros(N)
  368.     y[a] = c
  369.     for i in range(1,N):
  370.         y[i] = y[i-1] + (H * f(x[i-1],y[i-1]))
  371.     return x,y
  372.  
  373. x,y = EulerMethod(0,1,0,100,f) #approximate solution
  374.  
  375. plt.plot(x,y)
  376. plt.legend(["f(x)"])
  377. plt.title("Euler's Method of function approximation")
  378. plt.xlabel("x")
  379. plt.ylabel("f(x)")
  380.  
  381.  
  382. # 2. Test the function <b> EulerMethod </b> by using the data in Exercise 6
  383.  
  384. # In[ ]:
  385.  
  386.  
  387.  
  388.  
  389.  
  390. # # $$\text{Exercise 9}$$
  391. #
  392. # Consider a 25 L tank that will be filled with water in two different ways. In the first case, the water volume that enters the tank per time (rate of volume increase) is piecewise constant, while in the second case, it is continuously increasing.
  393. #
  394. # For each of these two cases, you are asked to develop a code that can compute how the total water volume V  in the tank will develop with time t over a period of 3 s. Your calculations must be based on the information given: the initial volume of water (1L in both cases), and the volume of water entering the tank per time. Write a python code to solve and plot the exact volume against the time for the rate $r(t)$ given by
  395. # 1. Case 1: piecewise constant
  396. # \begin{equation}
  397. # r(t) =
  398. #   \begin{cases}
  399. #     1 L/s      & \quad \text{if }\quad  0s<t<1s\\
  400. #     3 L/s  & \quad \text{if }\quad  1s\leq t< 2s\\
  401. #     7 L/s & \quad  \text{if }\quad  2s\leq t\leq 3s
  402. #   \end{cases}
  403. # \end{equation}
  404. #
  405. # 2. Case 2:  $r(t) = e^t$ for all $0< t\leq 3$.
  406. #
  407.  
  408. # In[131]:
  409.  
  410.  
  411. def r(t):
  412.         if 0<=t<1:
  413.             return t+1
  414.         elif 1<=t<2:
  415.             return 3*t+1
  416.         elif 2<=t<=3:
  417.             return 7*t+1
  418. vr = np.vectorize(r)
  419.  
  420. #def volume(a,b,t,N):
  421. #    v=np.zeros(N+1)
  422. #    for i in range(len(t)):
  423. #        return vr(i)
  424.  
  425.  
  426. a = 0
  427. b = 3
  428. N = 100
  429. t=np.linspace(a,b,N+1)
  430. v = vr(t)
  431. #v  = volume(a,b,t,N)
  432.  
  433.  
  434. # In[132]:
  435.  
  436.  
  437. def Vol_Exp(a,b,N):
  438.     t = np.linspace(a,b,N+1)
  439.     v_e = np.zeros(N+1)
  440.     v_e[0] = 1
  441.     for k in range(len(t)):
  442.         v_e[k]=np.exp(k)
  443.     return v_e
  444.  
  445. v_e = Vol_Exp(0,3,100)
  446. v_e
  447.  
  448.  
  449. # In[159]:
  450.  
  451.  
  452. def grapher(t,v,v_e,label):
  453.     plt.figure(figsize=(8,4))
  454.     plt.subplot(2,1,1)
  455.     plt.plot(t,v)
  456.     plt.title('Graph of Exact Values against time')
  457.     plt.xlabel('Time')
  458.     plt.ylabel('Exact Values')
  459.     plt.legend(['Piecewise'])
  460.     plt.show()
  461.     plt.figure(figsize=(12,6))
  462.     plt.subplot(2,2,1)
  463.     plt.plot(t,v_e)
  464.     plt.title('Graph of Exact Values against time')
  465.     plt.xlabel('Time')
  466.     plt.ylabel('Exact Values')
  467.     plt.legend(['Continuous'])
  468.     plt.show()
  469. grapher(t,v,v_e,['Graph of Exact Values against time','Exact Values','Piecewise','Continuous'])
  470.  
  471.  
  472. # In[ ]:
  473.  
  474.  
  475.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!