1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
|
\documentclass[12pt]{article}
\title{Solving Equations \& ODEs}
\author{FOSSEE}
\usepackage{listings}
\lstset{language=Python,
basicstyle=\ttfamily,
commentstyle=\itshape\bfseries,
showstringspaces=false,
}
\newcommand{\typ}[1]{\lstinline{#1}}
\usepackage[english]{babel}
\usepackage[latin1]{inputenc}
\usepackage{times}
\usepackage[T1]{fontenc}
\usepackage{ae,aecompl}
\usepackage{mathpazo,courier,euler}
\usepackage[scaled=.95]{helvet}
\usepackage{amsmath}
\begin{document}
\date{}
\vspace{-1in}
\begin{center}
\LARGE{Solving Equations \& ODEs}\\
\large{FOSSEE}
\end{center}
\section{Solving linear equations}
Consider following sets of equations:\\
\begin{align*}
3x + 2y - z & = 1 \\
2x - 2y + 4z & = -2 \\
-x + \frac{1}{2}y -z & = 0
\end{align*}\\
The matrix representation is:\\
\begin{center}
$A*x = B$
\end{center}
Where A is coefficient matrix(in this case 3x3)\\
B is constant matrix(1x3)\\
x is the required solution.\\
\begin{lstlisting}
In []: A = array([[3,2,-1], [2,-2,4], [-1, 0.5, -1]])
In []: B = array([[1], [-2], [0]])
In []: x = solve(A, B)
\end{lstlisting}
Solve the equation $A x = B$ for $x$.\\
To check whether solution is correct try this:
\begin{lstlisting}
In []: Ax = dot(A,x) #Matrix multiplication of A and x(LHS)
In []: allclose(Ax, B)
Out[]: True
\end{lstlisting}
\typ{allclose} Returns \typ{True} if two arrays(in above case Ax and B) are element-wise equal within a tolerance.
\newpage
\section{Finding roots}
\subsection{Polynomials}
\begin{center}
$x^2+6x+13=0$
\end{center}
to find roots, pylab provides \typ{roots} function.
\begin{lstlisting}
In []: coeffs = [1, 6, 13] #list of all coefficients
In []: roots(coeffs)
\end{lstlisting}
\subsection{functions}
Functions can be defined and used by following syntax:
\begin{lstlisting}
def func_name(arg1, arg2):
#function code
return ret_value
\end{lstlisting}
A simple example can be
\begin{lstlisting}
def expression(x):
y = x*sin(x)
return y
\end{lstlisting}
Above function when called with a argument, will return $xsin(x)$ value for that argument.
\begin{lstlisting}
In [95]: expression(pi/2)
Out[95]: 1.5707963267948966
In [96]: expression(pi/3)
Out[96]: 0.90689968211710881
\end{lstlisting}
\subsection{Roots of non-linear equations}
For Finding the roots of a non linear equation(defined as $f(x)=0$), around a starting estimate we use \typ{fsolve}:\\
\typ{In []: from scipy.optimize import fsolve}\\
\typ{fsolve} is not a part of \typ{pylab}, instead is a function in the \textbf{optimize} module of \textbf{scipy}, and hence we \textbf{import} it.\\
%\typ{fsolve} takes first argument as name of function, which evaluates $f(x)$, whose roots one wants to find. And second argument is starting estimate, around which roots are found.
For illustration, we want to find roots of equation:
\begin{center}
$f(x)=sin(x)+cos(x)^2$
\end{center}
So just like we did above, we define a function:
\begin{lstlisting}
In []: def f(x):
....: return sin(x)+cos(x)**2
....:
\end{lstlisting}
Now to find roots of this non linear equation, around a initial estimate value, say 0, we use \typ{fsolve} in following way:
\begin{lstlisting}
In []: fsolve(f, 0) #arguments are function name and estimate
Out[]: -0.66623943249251527
\end{lstlisting}
\section{ODE}
We wish to solve an (a system of) Ordinary Differential Equation. For this purpose, we shall use \typ{odeint}.\\
\typ{In []: from scipy.integrate import odeint}\\
\typ{odeint} is a function in the \textbf{integrate} module of \textbf{scipy}.\\
As an illustration, let us solve the ODE below:\\
$\frac{dy}{dt} = ky(L-y)$, L = 25000, k = 0.00003, y(0) = 250\\
We define a function (as below) that takes $y$ and time as arguments and returns the right hand side of the ODE.
\begin{lstlisting}
In []: def f(y, t):
.... k, L = 0.00003, 25000
.... return k*y*(L-y)
....
\end{lstlisting}
Next we define the time over which we wish to solve the ODE. We also note the initial conditions given to us.
\begin{lstlisting}
In []: t = linspace(0, 12, 61)
In []: y0 = 250
\end{lstlisting}
To solve the ODE, we call the \typ{odeint} function with three arguments - the function \typ{f}, initial conditions and the time vector.
\begin{lstlisting}
In []: y = odeint(f, y0, t)
\end{lstlisting}
Note: To solve a system of ODEs, we need to change the function to
return the right hand side of all the equations and the system and the
pass the required number of initial conditions to the \typ{odeint}
function.
\section{FFT}
\begin{itemize}
\item We have a simple signal $y(t)$
\item Find the FFT and plot it
\end{itemize}
\begin{lstlisting}
In []: t = linspace(0, 2*pi, 500)
In []: y = sin(4*pi*t)
In []: f = fft(y)
In []: freq = fftfreq(500, t[1] - t[0])
In []: plot(freq[:250], abs(f)[:250])
In []: grid()
\end{lstlisting}
\begin{itemize}
\item We now calculate the inverse Fourier transform.
\item Then, verify the solution obtained.
\end{itemize}
\begin{lstlisting}
In []: y1 = ifft(f) # inverse FFT
In []: allclose(y, y1)
Out[]: True
\end{lstlisting}
\begin{itemize}
\item Let us add some noise to the signal
\end{itemize}
\begin{lstlisting}
In []: yr = y + random(size=500)*0.2
In []: yn = y + normal(size=500)*0.2
In []: plot(t, yr)
In []: figure()
In []: plot(freq[:250],
...: abs(fft(yn))[:250])
\end{lstlisting}
\begin{itemize}
\item \typ{random}: produces uniform deviates in $[0, 1)$
\item \typ{normal}: draws random samples from a Gaussian
distribution
\item Useful to create a random matrix of any shape
\end{itemize}
\begin{itemize}
\item Now, we filter the noisy signal using a Wiener filter
\end{itemize}
\begin{lstlisting}
In []: from scipy import signal
In []: yc = signal.wiener(yn, 5)
In []: clf()
In []: plot(t, yc)
In []: figure()
In []: plot(freq[:250],
...: abs(fft(yc))[:250])
\end{lstlisting}
\section{Links and References}
\begin{itemize}
\item Documentation for Numpy and Scipy is available at:\\ http://docs.scipy.org/doc/
\item For "recipes" or worked examples of commonly-done tasks in SciPy explore: \\ http://www.scipy.org/Cookbook/
\end{itemize}
\end{document}
|
