\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}