122 \frametitle{Outline} |
122 \frametitle{Outline} |
123 \tableofcontents |
123 \tableofcontents |
124 % \pausesections |
124 % \pausesections |
125 \end{frame} |
125 \end{frame} |
126 |
126 |
127 \section{Integration} |
127 \section{Interpolation} |
128 |
128 |
129 \subsection{Quadrature} |
129 \begin{frame}[fragile] |
|
130 \frametitle{Interpolation} |
|
131 \begin{itemize} |
|
132 \item Let us begin with interpolation |
|
133 \item Let's use the L and T arrays and interpolate this data to obtain data at new points |
|
134 \end{itemize} |
|
135 \begin{lstlisting} |
|
136 In []: L = [] |
|
137 In []: T = [] |
|
138 In []: for line in open('pendulum.txt'): |
|
139 l, t = line.split() |
|
140 L.append(float(l)) |
|
141 T.append(float(t)) |
|
142 In []: L = array(L) |
|
143 In []: T = array(T) |
|
144 \end{lstlisting} |
|
145 \end{frame} |
|
146 |
|
147 %% \begin{frame}[fragile] |
|
148 %% \frametitle{Interpolation \ldots} |
|
149 %% \begin{small} |
|
150 %% \typ{In []: from scipy.interpolate import interp1d} |
|
151 %% \end{small} |
|
152 %% \begin{itemize} |
|
153 %% \item The \typ{interp1d} function returns a function |
|
154 %% \begin{lstlisting} |
|
155 %% In []: f = interp1d(L, T) |
|
156 %% \end{lstlisting} |
|
157 %% \item Functions can be assigned to variables |
|
158 %% \item This function interpolates between known data values to obtain unknown |
|
159 %% \end{itemize} |
|
160 %% \end{frame} |
|
161 |
|
162 %% \begin{frame}[fragile] |
|
163 %% \frametitle{Interpolation \ldots} |
|
164 %% \begin{lstlisting} |
|
165 %% In []: Ln = arange(0.1,0.99,0.005) |
|
166 %% # Interpolating! |
|
167 %% # The new values in range of old data |
|
168 %% In []: plot(L, T, 'o', Ln, f(Ln), '-') |
|
169 %% In []: f = interp1d(L, T, kind='cubic') |
|
170 %% # When kind not specified, it's linear |
|
171 %% # Others are ... |
|
172 %% # 'nearest', 'zero', |
|
173 %% # 'slinear', 'quadratic' |
|
174 %% \end{lstlisting} |
|
175 %% \end{frame} |
|
176 |
|
177 \begin{frame}[fragile] |
|
178 \frametitle{Spline Interpolation} |
|
179 \begin{small} |
|
180 \begin{lstlisting} |
|
181 In []: from scipy.interpolate import splrep |
|
182 In []: from scipy.interpolate import splev |
|
183 \end{lstlisting} |
|
184 \end{small} |
|
185 \begin{itemize} |
|
186 \item Involves two steps |
|
187 \begin{enumerate} |
|
188 \item Find out the spline curve, coefficients |
|
189 \item Evaluate the spline at new points |
|
190 \end{enumerate} |
|
191 \end{itemize} |
|
192 \end{frame} |
|
193 |
|
194 \begin{frame}[fragile] |
|
195 \frametitle{\typ{splrep}} |
|
196 To find the B-spline representation |
|
197 \begin{lstlisting} |
|
198 In []: tck = splrep(L, T) |
|
199 \end{lstlisting} |
|
200 Returns |
|
201 \begin{enumerate} |
|
202 \item the vector of knots, |
|
203 \item the B-spline coefficients |
|
204 \item the degree of the spline (default=3) |
|
205 \end{enumerate} |
|
206 \end{frame} |
|
207 |
|
208 \begin{frame}[fragile] |
|
209 \frametitle{\typ{splev}} |
|
210 To Evaluate a B-spline and it's derivatives |
|
211 \begin{lstlisting} |
|
212 In []: Lnew = arange(0.1,1,0.005) |
|
213 In []: Tnew = splev(Lnew, tck) |
|
214 |
|
215 #To obtain derivatives of the spline |
|
216 #use der=1, 2,.. for 1st, 2nd,.. order |
|
217 In []: Tnew = splev(Lnew, tck, der=1) |
|
218 \end{lstlisting} |
|
219 \end{frame} |
|
220 |
|
221 %% \begin{frame}[fragile] |
|
222 %% \frametitle{Interpolation \ldots} |
|
223 %% \begin{itemize} |
|
224 %% \item |
|
225 %% \end{itemize} |
|
226 %% \end{frame} |
|
227 |
|
228 \section{Differentiation} |
|
229 |
|
230 \begin{frame}[fragile] |
|
231 \frametitle{Numerical Differentiation} |
|
232 \begin{itemize} |
|
233 \item Given function $f(x)$ or data points $y=f(x)$ |
|
234 \item We wish to calculate $f^{'}(x)$ at points $x$ |
|
235 \item Taylor series - finite difference approximations |
|
236 \end{itemize} |
|
237 \begin{center} |
|
238 \begin{tabular}{l l} |
|
239 $f(x+h)=f(x)+h.f^{'}(x)$ &Forward \\ |
|
240 $f(x-h)=f(x)-h.f^{'}(x)$ &Backward |
|
241 \end{tabular} |
|
242 \end{center} |
|
243 \end{frame} |
|
244 |
|
245 \begin{frame}[fragile] |
|
246 \frametitle{Forward Difference} |
|
247 \begin{lstlisting} |
|
248 In []: x = linspace(0, 2*pi, 100) |
|
249 In []: y = sin(x) |
|
250 In []: deltax = x[1] - x[0] |
|
251 \end{lstlisting} |
|
252 Obtain the finite forward difference of y |
|
253 \end{frame} |
|
254 |
|
255 \begin{frame}[fragile] |
|
256 \frametitle{Forward Difference \ldots} |
|
257 \begin{lstlisting} |
|
258 In []: fD = (y[1:] - y[:-1]) / deltax |
|
259 In []: plot(x, y, x[:-1], fD) |
|
260 \end{lstlisting} |
|
261 \begin{center} |
|
262 \includegraphics[height=2in, interpolate=true]{data/fwdDiff} |
|
263 \end{center} |
|
264 \end{frame} |
|
265 |
|
266 \begin{frame}[fragile] |
|
267 \frametitle{Example} |
|
268 \begin{itemize} |
|
269 \item Given x, y positions of a particle in \typ{pos.txt} |
|
270 \item Find velocity \& acceleration in x, y directions |
|
271 \end{itemize} |
|
272 \small{ |
|
273 \begin{center} |
|
274 \begin{tabular}{| c | c | c |} |
|
275 \hline |
|
276 $X$ & $Y$ \\ \hline |
|
277 0. & 0.\\ \hline |
|
278 0.25 & 0.47775\\ \hline |
|
279 0.5 & 0.931\\ \hline |
|
280 0.75 & 1.35975\\ \hline |
|
281 1. & 1.764\\ \hline |
|
282 1.25 & 2.14375\\ \hline |
|
283 \vdots & \vdots\\ \hline |
|
284 \end{tabular} |
|
285 \end{center}} |
|
286 \end{frame} |
|
287 |
|
288 \begin{frame}[fragile] |
|
289 \frametitle{Example \ldots} |
|
290 \begin{itemize} |
|
291 \item Read the file |
|
292 \item Obtain an array of x, y |
|
293 \item Obtain velocity and acceleration |
|
294 \item use \typ{deltaT = 0.05} |
|
295 \end{itemize} |
|
296 \begin{lstlisting} |
|
297 In []: X = [] |
|
298 In []: Y = [] |
|
299 In []: for line in open('location.txt'): |
|
300 .... points = line.split() |
|
301 .... X.append(float(points[0])) |
|
302 .... Y.append(float(points[1])) |
|
303 In []: S = array([X, Y]) |
|
304 \end{lstlisting} |
|
305 \end{frame} |
|
306 |
|
307 |
|
308 \begin{frame}[fragile] |
|
309 \frametitle{Example \ldots} |
|
310 \begin{itemize} |
|
311 \item use \typ{deltaT = 0.05} |
|
312 \end{itemize} |
|
313 \begin{lstlisting} |
|
314 In []: deltaT = 0.05 |
|
315 |
|
316 In []: v = (S[:,1:]-S[:,:-1])/deltaT |
|
317 |
|
318 In []: a = (v[:,1:]-v[:,:-1])/deltaT |
|
319 \end{lstlisting} |
|
320 Try Plotting the position, velocity \& acceleration. |
|
321 \end{frame} |
|
322 |
|
323 \section{Quadrature} |
130 |
324 |
131 \begin{frame}[fragile] |
325 \begin{frame}[fragile] |
132 \frametitle{Quadrature} |
326 \frametitle{Quadrature} |
133 \begin{itemize} |
327 \begin{itemize} |
134 \item We wish to find area under a curve |
328 \item We wish to find area under a curve |
135 \item Area under $(sin(x) + x^2)$ in $(0,1)$ |
329 \item Area under $(sin(x) + x^2)$ in $(0,1)$ |
136 \item scipy has functions to do that |
330 \item scipy has functions to do that |
137 \end{itemize} |
331 \end{itemize} |
138 \small{\typ{In []: from scipy.integrate import quad}} |
332 \begin{small} |
|
333 \typ{In []: from scipy.integrate import quad} |
|
334 \end{small} |
139 \begin{itemize} |
335 \begin{itemize} |
140 \item Inputs - function to integrate, limits |
336 \item Inputs - function to integrate, limits |
141 \end{itemize} |
337 \end{itemize} |
142 \begin{lstlisting} |
338 \begin{lstlisting} |
143 In []: x = 0 |
339 In []: x = 0 |
144 In []: quad(sin(x)+x**2, 0, 1) |
340 In []: quad(sin(x)+x**2, 0, 1) |
145 \end{lstlisting} |
341 \end{lstlisting} |
|
342 \begin{small} |
146 \alert{\typ{error:}} |
343 \alert{\typ{error:}} |
147 \typ{First argument must be a callable function.} |
344 \typ{First argument must be a callable function.} |
|
345 \end{small} |
148 \end{frame} |
346 \end{frame} |
149 |
347 |
150 \begin{frame}[fragile] |
348 \begin{frame}[fragile] |
151 \frametitle{Functions - Definition} |
349 \frametitle{Functions - Definition} |
|
350 We have been using them all along. Now let's see how to define them. |
152 \begin{lstlisting} |
351 \begin{lstlisting} |
153 In []: def f(x): |
352 In []: def f(x): |
154 return sin(x)+x**2 |
353 return sin(x)+x**2 |
155 In []: quad(f, 0, 1) |
354 In []: quad(f, 0, 1) |
156 \end{lstlisting} |
355 \end{lstlisting} |
157 \begin{itemize} |
356 \begin{itemize} |
158 \item \typ{def} |
357 \item \typ{def} |
|
358 \item name |
159 \item arguments |
359 \item arguments |
160 \item \typ{return} |
360 \item \typ{return} |
161 \end{itemize} |
361 \end{itemize} |
162 \end{frame} |
362 \end{frame} |
163 |
363 |
173 In []: f(0) |
373 In []: f(0) |
174 Out[]: 0.0 |
374 Out[]: 0.0 |
175 In []: f(1) |
375 In []: f(1) |
176 Out[]: 1.8414709848078965 |
376 Out[]: 1.8414709848078965 |
177 \end{lstlisting} |
377 \end{lstlisting} |
178 \end{frame} |
378 More on Functions later \ldots |
179 |
|
180 |
|
181 \begin{frame}[fragile] |
|
182 \frametitle{Functions - Default Arguments} |
|
183 \begin{lstlisting} |
|
184 In []: def f(x=1): |
|
185 return sin(x)+x**2 |
|
186 In []: f(10) |
|
187 Out[]: 99.455978889110625 |
|
188 In []: f(1) |
|
189 Out[]: 1.8414709848078965 |
|
190 In []: f() |
|
191 Out[]: 1.8414709848078965 |
|
192 \end{lstlisting} |
|
193 \end{frame} |
|
194 |
|
195 \begin{frame}[fragile] |
|
196 \frametitle{Functions - Keyword Arguments} |
|
197 \begin{lstlisting} |
|
198 In []: def f(x=1, y=pi): |
|
199 return sin(y)+x**2 |
|
200 In []: f() |
|
201 Out[]: 1.0000000000000002 |
|
202 In []: f(2) |
|
203 Out[]: 4.0 |
|
204 In []: f(y=2) |
|
205 Out[]: 1.9092974268256817 |
|
206 In []: f(y=pi/2,x=0) |
|
207 Out[]: 1.0 |
|
208 \end{lstlisting} |
|
209 \end{frame} |
|
210 |
|
211 \begin{frame}[fragile] |
|
212 \frametitle{More on functions} |
|
213 \begin{itemize} |
|
214 \item Scope of variables in the function is local |
|
215 \item Mutable items are \alert{passed by reference} |
|
216 \item First line after definition may be a documentation string |
|
217 (\alert{recommended!}) |
|
218 \item Function definition and execution defines a name bound to the |
|
219 function |
|
220 \item You \emph{can} assign a variable to a function! |
|
221 \end{itemize} |
|
222 \end{frame} |
379 \end{frame} |
223 |
380 |
224 \begin{frame}[fragile] |
381 \begin{frame}[fragile] |
225 \frametitle{Quadrature \ldots} |
382 \frametitle{Quadrature \ldots} |
226 \begin{lstlisting} |
383 \begin{lstlisting} |
227 In []: quad(f, 0, 1) |
384 In []: quad(f, 0, 1) |
228 \end{lstlisting} |
385 \end{lstlisting} |
229 Returns the integral and an estimate of the absolute error in the result. |
386 Returns the integral and an estimate of the absolute error in the result. |
230 \begin{itemize} |
387 \begin{itemize} |
231 \item Use \typ{dblquad} for Double integrals |
388 \item Look at \typ{dblquad} for Double integrals |
232 \item Use \typ{tplquad} for Triple integrals |
389 \item Use \typ{tplquad} for Triple integrals |
233 \end{itemize} |
390 \end{itemize} |
234 \end{frame} |
|
235 |
|
236 \subsection{ODEs} |
|
237 |
|
238 \begin{frame}[fragile] |
|
239 \frametitle{ODE Integration} |
|
240 We shall use the simple ODE of a simple pendulum. |
|
241 \begin{equation*} |
|
242 \ddot{\theta} = -\frac{g}{L}sin(\theta) |
|
243 \end{equation*} |
|
244 \begin{itemize} |
|
245 \item This equation can be written as a system of two first order ODEs |
|
246 \end{itemize} |
|
247 \begin{align} |
|
248 \dot{\theta} &= \omega \\ |
|
249 \dot{\omega} &= -\frac{g}{L}sin(\theta) \\ |
|
250 \text{At}\ t &= 0 : \nonumber \\ |
|
251 \theta = \theta_0\quad & \&\quad \omega = 0 \nonumber |
|
252 \end{align} |
|
253 \end{frame} |
|
254 |
|
255 \begin{frame}[fragile] |
|
256 \frametitle{Solving ODEs using SciPy} |
|
257 \begin{itemize} |
|
258 \item We use the \typ{odeint} function from scipy to do the integration |
|
259 \item Define a function as below |
|
260 \end{itemize} |
|
261 \begin{lstlisting} |
|
262 In []: def pend_int(unknown, t, p): |
|
263 .... theta, omega = unknown |
|
264 .... g, L = p |
|
265 .... f=[omega, -(g/L)*sin(theta)] |
|
266 .... return f |
|
267 .... |
|
268 \end{lstlisting} |
|
269 \end{frame} |
|
270 |
|
271 \begin{frame}[fragile] |
|
272 \frametitle{Solving ODEs using SciPy \ldots} |
|
273 \begin{itemize} |
|
274 \item \typ{t} is the time variable \\ |
|
275 \item \typ{p} has the constants \\ |
|
276 \item \typ{initial} has the initial values |
|
277 \end{itemize} |
|
278 \begin{lstlisting} |
|
279 In []: t = linspace(0, 10, 101) |
|
280 In []: p=(-9.81, 0.2) |
|
281 In []: initial = [10*2*pi/360, 0] |
|
282 \end{lstlisting} |
|
283 \end{frame} |
|
284 |
|
285 \begin{frame}[fragile] |
|
286 \frametitle{Solving ODEs using SciPy \ldots} |
|
287 |
|
288 \small{\typ{In []: from scipy.integrate import odeint}} |
|
289 \begin{lstlisting} |
|
290 In []: pend_sol = odeint(pend_int, |
|
291 initial,t, |
|
292 args=(p,)) |
|
293 \end{lstlisting} |
|
294 \end{frame} |
391 \end{frame} |
295 |
392 |
296 \begin{frame} |
393 \begin{frame} |
297 \frametitle{Things we have learned} |
394 \frametitle{Things we have learned} |
298 \begin{itemize} |
395 \begin{itemize} |
|
396 \item Interpolation |
|
397 \item Differentiation |
299 \item Functions |
398 \item Functions |
300 \begin{itemize} |
399 \begin{itemize} |
301 \item Definition |
400 \item Definition |
302 \item Calling |
401 \item Calling |
303 \item Default Arguments |
402 \item Default Arguments |