Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #!/usr/bin/env python
- # coding: utf-8
- # In[31]:
- get_ipython().run_cell_magic('javascript', '', 'MathJax.Hub.Config({\n TeX: { equationNumbers: { autoNumber: "AMS" } }\n});')
- # # $$\text{Exercise 1}$$
- # 1. Let us consider the array $M$ defined below
- # In[47]:
- import numpy as np
- import matplotlib.pyplot as plt
- M = np.array([[ 80842, 333008, 202553, 140037, 81969, 63857, 42105, 261540],
- [481981, 176739, 489984, 326386, 110795, 394863, 25024, 38317],
- [ 49982, 408830, 485118, 16119, 407675, 231729, 265455, 109413],
- [103399, 174677, 343356, 301717, 224120, 401101, 140473, 254634],
- [112262, 25063, 108262, 375059, 406983, 208947, 115641, 296685],
- [444899, 129585, 171318, 313094, 425041, 188411, 335140, 141681],
- [ 59641, 211420, 287650, 8973, 477425, 382803, 465168, 3975],
- [ 32213, 160603, 275485, 388234, 246225, 56174, 244097, 9350],
- [496966, 225516, 273338, 73335, 283013, 212813, 38175, 282399],
- [318413, 337639, 379802, 198049, 101115, 419547, 260219, 325793],
- [148593, 425024, 348570, 117968, 107007, 52547, 180346, 178760],
- [305186, 262153, 11835, 449971, 494184, 472031, 353049, 476442],
- [ 35455, 191553, 384154, 29917, 187599, 68912, 428968, 69638],
- [ 20140, 220691, 163112, 47862, 474879, 15713, 279173, 161443],
- [ 76851, 450401, 193338, 189948, 211826, 384657, 17636, 189954],
- [182580, 403458, 63004, 87202, 17143, 308769, 316192, 408259],
- [379223, 322891, 302814, 309594, 13650, 68095, 327748, 27138],
- [ 33346, 202960, 221234, 61986, 96524, 230606, 266487, 299766],
- [ 46264, 33510, 434647, 26750, 412754, 346058, 263454, 270444],
- [182770, 486269, 402948, 105055, 463911, 378799, 329053, 16793],
- [379586, 293162, 121602, 144605, 105553, 71148, 40126, 183534],
- [493338, 5241, 498084, 189451, 219466, 12869, 482845, 16529],
- [145955, 461866, 427009, 70841, 318111, 199631, 99330, 407887],
- [133826, 214053, 295543, 339707, 199545, 194494, 304822, 337793],
- [395524, 206326, 308488, 297633, 472409, 487606, 422244, 443487],
- [464749, 75275, 299486, 248984, 76076, 197544, 75103, 394086],
- [404365, 67309, 267556, 96483, 166099, 203867, 57120, 211249],
- [ 86992, 286426, 130451, 419955, 109839, 332663, 101095, 50820],
- [108699, 436141, 333670, 45890, 425798, 6967, 155851, 360687],
- [252907, 121138, 81509, 367876, 471831, 226950, 349068, 197092]])
- M
- # 1. What is the average value of the second and fourth columns?
- # In[38]:
- np.average(M[:,1:4:2])
- # In[8]:
- np.average(M[:,1:4:2])
- # 2. What is the average value of the last 5 rows of the second and fourth columns?
- # In[10]:
- np.average(M[-5:,1:4:2])
- # 3. Update $M$ by Replacing all its even integers by 0.
- # In[11]:
- M[M%2==0]=0
- M
- # # $$\text{Exercise 2}$$
- #
- # 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}$),
- # then we have
- # $$\int_a^bf(x)dx \approx \sum_{k=1}^N \frac{f(x_{k-1})+f(x_k)}{2}H = A$$
- #
- #
- # Write a function named <b>Trap</b> that takes $a,b,N, f$ as inputs and return A
- # In[39]:
- import numpy as np
- def Trap(a,b,N,f):
- x = np.linspace(a,b,N+1)
- H = x[1]-x[0]
- A = 0
- for k in range (1,N+1):
- A+=(H/2)*(f(x[k-1])+f(x[k]))
- return A
- Trap(0,10,20,lambda x:x**2)
- # # $$\text{Exercise 3}$$
- #
- # Let us consider the function $f$ defined on the set of integer as follows
- #
- # \begin{equation}
- # f(n) =
- # \begin{cases}
- # n/2 & \quad \text{if } n \text{ is even}\\
- # -(n)/3 +1 & \quad \text{if } n \text{ is odd and divisible by 3}\\
- # (n+1)/2 +1 & \quad \text{else}
- # \end{cases}
- # \end{equation}
- #
- # 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]$.
- # In[40]:
- def f(n):
- if n%2==0:
- return n/2
- elif n%2==1 and n%3==0:
- return (-(n)/3)+1
- else:
- return ((n+1)/2)+1
- vf = np.vectorize(f)
- def CompFunc(a,b):
- if a<b:
- nlist = np.arange(a,b+1)
- flist = vf(nlist)
- return flist
- CompFunc(0,9)
- # # $$\text{Exercise 4}$$
- #
- # Write a python code to create the following numpy arrays
- # $$
- # A = \begin{pmatrix}
- # 1 & 4 & 6 \\
- # 0 & -3 & 2 \\
- # -2 & -2 & -2
- # \end{pmatrix}, \quad
- # B = \begin{pmatrix}
- # 2 & -1 & 0 \\
- # 2 & -1 & 0 \\
- # 2 & -3 & 1
- # \end{pmatrix},
- # $$
- # and compute
- # - $A-B$,
- # - $4A + B$,
- # - $trace(A)$,
- # - $B^t$ the transpose of $B$,
- # - $AB$,
- # - $BA^t$,
- # - the determinant of $A$
- # In[13]:
- A = np.array([[1,4,6],[0,-3,2],[-2,-2,-2]])
- B = np.array([[2,-1,0],[2,-1,0],[2,-3,1]])
- # In[14]:
- A-B
- # In[15]:
- 4*A+B
- # In[16]:
- np.trace(A)
- # In[17]:
- B.T
- # In[ ]:
- A
- # In[ ]:
- # # $$\text{Exercise 5}$$
- #
- # write a python code to solve the system of equation
- # \begin{equation}\label{sysEqu}
- # \left\lbrace
- # \begin{array}{ll}
- # 2x + y -2z &= 3,\\
- # x-y-z &= 0,\\
- # x+y+3z &=12
- # \end{array}\right.
- # \end{equation}
- # In[42]:
- A = np.array([[2,1,-2],[1,-1,-1],[1,1,3]])
- B = np.array([3,0,12])
- P = np.linalg.solve(A,B)
- P
- # # $$\text{Exercise 6}$$
- #
- # Let us consider the sequence $U_n$ given by
- # \begin{equation}\label{fib}
- # \left\lbrace
- # \begin{array}{ll}
- # U_0 &= 1,\\
- # U_1 &= 2,\\
- # U_{n} &=-3U_{n-1} +U_{n-2}, \;\; \forall\; n=2,3,4\cdots
- # \end{array}\right.
- # \end{equation}
- #
- # 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}.
- # In[118]:
- def SeqTerms(n):
- U = np.array([])
- if(n >= 0):
- if(n == 0):
- return U
- elif(n == 1):
- return np.append(U,1)
- elif(n == 2):
- return np.append(U,[1,2])
- else:
- U = np.append(U,[1,2])
- for i in range(2,n):
- U = np.append(U,-3*U[i-1] + U[i-2])
- return U.astype(int)
- SeqTerms(10)
- # In[71]:
- len([])
- # # $$\text{Exercise 7}$$
- #
- # 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}$).
- #
- # 1. Write a python code that use Euler’s method to solve the initial value problem
- # \begin{equation}\label{eul1}
- # \begin{cases}
- # y' = \dfrac{2-2xy}{x^2+1}, & \quad \text{on } [0, 1]\\\\
- # y(0) = 1,
- # \end{cases}
- # \end{equation}
- # i.e. generates array of $y_k\approx g(x_k)=g_k$ where $g$ is the exact solution of \eqref{eul1}
- #
- # 2. The exact solution of the initial value problem is given by
- # \begin{equation}\label{exact}
- # g(x) = \dfrac{2x+1}{x^2+1}
- # \end{equation}
- #
- # use $subplot$ to plot side by side
- # - the exact and appromixate solution in the same window (i.e. ($y_k$ vs $x_k$ and $g_k$ vs $x_k$))
- # - the error $e_k = |y_k -g_k| $ against $x_k$.
- # In[108]:
- f = lambda x,y: (2-(2*x*y))/(x**2 +1)
- g = lambda x: (2*x +1)/(x**2 +1)
- z = lambda f,g: abs(f-g)
- def Euler(a,b,n,f):
- x = np.linspace(a,b,n)
- H = (b-a)/n
- y = np.zeros(n)
- y[0] = 1
- for i in range(1,n):
- y[i] = y[i-1] + (H * f(x[i-1],y[i-1]))
- return x,y
- # In[135]:
- x,y = Euler(0,1,100,f) #approximate solution
- plt.plot(x,y)
- y_0 = g(x) #actual solution
- #plt.plot(x,y_0)
- #plt.legend(["f(x)","g(x)"])
- #plt.title("Euler's Method of function approximation")
- #plt.xlabel("x")
- #plt.ylabel("f(x),g(x)")
- err = z(y,y_0)
- #print(err)
- fig, (ax1, ax2) = plt.subplots(1, 2)
- fig.suptitle('Function approximation')
- ax1.plot(x, y)
- ax1.title.set_text('Euler\'s Method of function approximation')
- ax2.plot(x, err)
- ax2.title.set_text('Approximation Error')
- # In[70]:
- x = np.linspace(0,10,100)
- y = np.exp(x)
- z = np.exp(-x) + np.sin(x)
- fig, (ax1, ax2) = plt.subplots(1, 2)
- fig.suptitle('Horizontally stacked subplots')
- ax1.plot(x, y)
- ax1.title.set_text('First Plot')
- ax2.plot(x, z)
- ax2.title.set_text('Second Plot')
- # # $$\text{Exercise 8}\quad (\text{generalization of previous exercise})$$
- #
- # 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
- #
- # \begin{equation}\label{eul2}
- # \begin{cases}
- # y' = f(x,y), & \quad \text{on } [a, b]\\\\
- # y(a) = c,
- # \end{cases}
- # \end{equation}
- #
- # 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.
- # In[116]:
- def EulerMethod(a,b,c,N,f):
- x = np.linspace(a,b,N)
- H = (b-a)/N
- y = np.zeros(N)
- y[a] = c
- for i in range(1,N):
- y[i] = y[i-1] + (H * f(x[i-1],y[i-1]))
- return x,y
- x,y = EulerMethod(0,1,0,100,f) #approximate solution
- plt.plot(x,y)
- plt.legend(["f(x)"])
- plt.title("Euler's Method of function approximation")
- plt.xlabel("x")
- plt.ylabel("f(x)")
- # 2. Test the function <b> EulerMethod </b> by using the data in Exercise 6
- # In[ ]:
- # # $$\text{Exercise 9}$$
- #
- # 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.
- #
- # 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
- # 1. Case 1: piecewise constant
- # \begin{equation}
- # r(t) =
- # \begin{cases}
- # 1 L/s & \quad \text{if }\quad 0s<t<1s\\
- # 3 L/s & \quad \text{if }\quad 1s\leq t< 2s\\
- # 7 L/s & \quad \text{if }\quad 2s\leq t\leq 3s
- # \end{cases}
- # \end{equation}
- #
- # 2. Case 2: $r(t) = e^t$ for all $0< t\leq 3$.
- #
- # In[131]:
- def r(t):
- if 0<=t<1:
- return t+1
- elif 1<=t<2:
- return 3*t+1
- elif 2<=t<=3:
- return 7*t+1
- vr = np.vectorize(r)
- #def volume(a,b,t,N):
- # v=np.zeros(N+1)
- # for i in range(len(t)):
- # return vr(i)
- a = 0
- b = 3
- N = 100
- t=np.linspace(a,b,N+1)
- v = vr(t)
- #v = volume(a,b,t,N)
- # In[132]:
- def Vol_Exp(a,b,N):
- t = np.linspace(a,b,N+1)
- v_e = np.zeros(N+1)
- v_e[0] = 1
- for k in range(len(t)):
- v_e[k]=np.exp(k)
- return v_e
- v_e = Vol_Exp(0,3,100)
- v_e
- # In[159]:
- def grapher(t,v,v_e,label):
- plt.figure(figsize=(8,4))
- plt.subplot(2,1,1)
- plt.plot(t,v)
- plt.title('Graph of Exact Values against time')
- plt.xlabel('Time')
- plt.ylabel('Exact Values')
- plt.legend(['Piecewise'])
- plt.show()
- plt.figure(figsize=(12,6))
- plt.subplot(2,2,1)
- plt.plot(t,v_e)
- plt.title('Graph of Exact Values against time')
- plt.xlabel('Time')
- plt.ylabel('Exact Values')
- plt.legend(['Continuous'])
- plt.show()
- grapher(t,v,v_e,['Graph of Exact Values against time','Exact Values','Piecewise','Continuous'])
- # In[ ]:
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement