Difference between revisions of "Numerical Methods for Ordinary Differential Equations"

From AstroBaki
Jump to navigationJump to search
Line 70: Line 70:
 
\\[12pt]
 
\\[12pt]
 
\\[12pt]
 
\\[12pt]
 +
 +
$\textbf{Example \, 1:}$
  
 
Now that we have a simple integration scheme, let's see how we'd execute this in python. Say our ODE is:
 
Now that we have a simple integration scheme, let's see how we'd execute this in python. Say our ODE is:
Line 79: Line 81:
 
\begin{verbatim}
 
\begin{verbatim}
  
.import numpy as np
+
. import numpy as np
.import matplotlib.pyplot as plt
+
. import matplotlib.pyplot as plt
 
.
 
.
.# plot analytic solution
+
. # plot analytic solution
.plt.figure(figsize=(10,4))
+
. plt.figure(figsize=(10,4))
.x_analytic = np.linspace(0,5,100)     
+
. x_analytic = np.linspace(0,5,100)     
.y_analytic = np.exp(x_analytic)
+
. y_analytic = np.exp(x_analytic)
.plt.plot(x_analytic, y_analytic, 'k--', label='analytic solution')
+
. plt.plot(x_analytic, y_analytic, 'k--', label='analytic solution')
 
.
 
.
.# define initial condition and values of h that will be used  
+
. # define initial condition and values of h that will be used  
.x0 = 0
+
. x0 = 0
.y0 = 1
+
. y0 = 1
.h_values = np.array([0.1,0.5,1.0])
+
. h_values = np.array([0.1,0.5,1.0])
 
.
 
.
.for h in h_values:
+
. for h in h_values:
.   # define arrays that can be used to record x and y
+
.     # define arrays that can be used to record x and y
.   x = np.zeros([100])
+
.     x = np.zeros([100])
.   x[0] = x0
+
.     x[0] = x0
.   y = np.zeros([100])
+
.     y = np.zeros([100])
.   y[0] = y0
+
.     y[0] = y0
.   # execute integration for current value of h
+
.     # execute integration for current value of h
.   for i in range(1, len(y)):
+
.     for i in range(1, len(y)):
.           x[i] = x[i-1] + h
+
.             x[i] = x[i-1] + h
.           dydx = y[i-1]
+
.             dydx = y[i-1]
.           y[i] = y[i-1] + h*dydx
+
.             y[i] = y[i-1] + h*dydx
.   # plot curve for current value of h
+
.     # plot curve for current value of h
.   plt.plot(np.trim_zeros(x, 'b'), np.trim_zeros(y, 'b'), label='h = '+str(h))
+
.     plt.plot(np.trim_zeros(x, 'b'), np.trim_zeros(y, 'b'), label='h = '+str(h))
 
.
 
.
.# remainder of plotting routine
+
. # remainder of plotting routine
.plt.xlabel('x', fontsize=12)
+
. plt.xlabel('x', fontsize=12)
.plt.ylabel('y', fontsize=12)
+
. plt.ylabel('y', fontsize=12)
.plt.xlim([0,7])
+
. plt.xlim([0,7])
.plt.ylim([0,100])
+
. plt.ylim([0,100])
.plt.legend(loc=2, fontsize=12)
+
. plt.legend(loc=2, fontsize=12)
.plt.tight_layout()
+
. plt.tight_layout()
.plt.show()
+
. plt.show()
  
 
\end{verbatim}
 
\end{verbatim}
Line 127: Line 129:
 
\\[12pt]
 
\\[12pt]
  
Since this is an astrophysics course, let's conclude this section with a more relevant example. Let's calculate the trajectory of a relativistic electron in the presence of an electric field (see: Synchrotron radiation). In this case,  
+
$\textbf{Example \, 2:}$
 +
 
 +
Since this is an astrophysics course, we will conclude this section with a more relevant example. Let's calculate the trajectory of a relativistic electron in the presence of a magnetic field (see: [https://casper.ssl.berkeley.edu/astrobaki/index.php/Synchrotron_Radiation Synchrotron Radiation]). Assume the electron begins at the origin and has an initial velocity of $\vec{v} = (0.99c, 0, 0.99c)$. Also assume a magnetic field with a strength of 1e-4 Gauss that points along the x-axis, or $\vec{B} = (\mathrm{1e-4}, 0, 0) \, \mathrm{Gauss}$. In this case, we need to apply Euler's method twice within each iteration of the integration. First, we calculate the new position of the electron given its instantaneous velocity:
 +
 
 +
$$\vec{x}_{i+1} = \vec{x}_{i} + h \vec{v}_{i}$$
 +
 
 +
and then we calculate the new velocity of the electron given its instantaneous acceleration:
 +
 
 +
$$\vec{v}_{i+1} = \vec{v}_{i} + h \vec{a}_{i} = \vec{v}_{i} + h \frac{e}{m_e} \left( \frac{\vec{v}_{i}}{c} \times \vec{B} \right)$$.
 +
 
 +
In our code, we will use a step size of 1e-6 seconds and integrate for a total of 0.1 second.
  
 
\subsection*{ The Runge-Kutta Method }
 
\subsection*{ The Runge-Kutta Method }

Revision as of 17:53, 11 December 2018

Short Topical Videos

Reference Materials


<latex> \documentclass[11pt]{article} \def\inv#1Template:1 \over \def\ddtTemplate:D \over dt \def\mean#1{\left\langle {#1}\right\rangle} \def\sigot{\sigma_{12}} \def\sigto{\sigma_{21}} \def\eval#1{\big|_{#1}} \def\tr{\nabla} \def\dce{\vec\tr\times\vec E} \def\dcb{\vec\tr\times\vec B} \def\wz{\omega_0} \def\ef{\vec E} \def\ato{{A_{21}}} \def\bto{{B_{21}}} \def\bot{{B_{12}}} \def\bfieldTemplate:\vec B \def\apTemplate:A^\prime \def\xp{{x^{\prime}}} \def\yp{{y^{\prime}}} \def\zp{{z^{\prime}}} \def\tp{{t^{\prime}}} \def\upxTemplate:U x^\prime \def\upyTemplate:U y^\prime \def\e#1{\cdot10^{#1}} \def\hf{\frac12} \def\^{\hat } \def\.{\dot } \usepackage{fullpage} \usepackage{amsmath} \usepackage{eufrak}

\begin{document} \section{ Numerical Integration Methods }

Numerical integration is used to solve ODEs that cannot be solved analytically, generally through discretization of the ODE. Since the conception of the modern computer, numerical integration methods have become an essential tool in the physical sciences and beyond. Here we will describe two of many such methods, along with some sample code.

\subsection*{ The Euler Method }

The Euler Method is the simplest numerical integration method, and is most likely the only one you will need to use in C207. Here, we'll walk through how it works and show some examples of how to implement it. Imagine you have some ODE:

$$\frac{dy}{dx} = f(x,y)$$

with some initial condition $y(x_0) = y_0$. Now, let's say we want to determine the value of $y_1$ at some new value of $x_1$ that is close to $x_0$. We can estimate this by rearranging the ODE above into the following form:

$$\frac{y_{1} - y_{0}}{x_{1} - x_{0}} = f(x_0,y_0)$$

which we can write as:

$$y_{1} = y_{0} + (x_{1} - x_{0}) f(x_0,y_0) = y_{0} + h f(x_0,y_0)$$

where $(x_{1} - x_{0}) = h$ is called the "step size. We can generalize this equation with the following recursive relation:

$$ y_{i+1} = y_{i} + h f(x_{i},y_{i}) $$

where the step size is typically set as some constant value. As we will see, smaller step sizes produce more accurate integration results. More precisely, the error is on the order of $O(h)$ (i.e., reducing $h$ by a factor of two makes the Euler method twice as accurate).

\\[12pt] \\[12pt]

$\textbf{Example \, 1:}$

Now that we have a simple integration scheme, let's see how we'd execute this in python. Say our ODE is:

$$\frac{dy}{dx} = y$$

and our initial condition is $y(0) = 1$. We can solve this analytically, but let's see what we get using the Euler method for different values of $h$.

\begin{verbatim}

. import numpy as np . import matplotlib.pyplot as plt . . # plot analytic solution . plt.figure(figsize=(10,4)) . x_analytic = np.linspace(0,5,100) . y_analytic = np.exp(x_analytic) . plt.plot(x_analytic, y_analytic, 'k--', label='analytic solution') . . # define initial condition and values of h that will be used . x0 = 0 . y0 = 1 . h_values = np.array([0.1,0.5,1.0]) . . for h in h_values: . # define arrays that can be used to record x and y . x = np.zeros([100]) . x[0] = x0 . y = np.zeros([100]) . y[0] = y0 . # execute integration for current value of h . for i in range(1, len(y)): . x[i] = x[i-1] + h . dydx = y[i-1] . y[i] = y[i-1] + h*dydx . # plot curve for current value of h . plt.plot(np.trim_zeros(x, 'b'), np.trim_zeros(y, 'b'), label='h = '+str(h)) . . # remainder of plotting routine . plt.xlabel('x', fontsize=12) . plt.ylabel('y', fontsize=12) . plt.xlim([0,7]) . plt.ylim([0,100]) . plt.legend(loc=2, fontsize=12) . plt.tight_layout() . plt.show()

\end{verbatim}

\begin{center} Numerical Euler1.png \end{center}

As we can see, the curve obtained using the Euler method converges to the analytic solution as $h \rightarrow 0$. Therefore, you must be careful to use a value of $h$ small enough to produce realistic results when performing numerical integrations.

\\[12pt] \\[12pt]

$\textbf{Example \, 2:}$

Since this is an astrophysics course, we will conclude this section with a more relevant example. Let's calculate the trajectory of a relativistic electron in the presence of a magnetic field (see: Synchrotron Radiation). Assume the electron begins at the origin and has an initial velocity of $\vec{v} = (0.99c, 0, 0.99c)$. Also assume a magnetic field with a strength of 1e-4 Gauss that points along the x-axis, or $\vec{B} = (\mathrm{1e-4}, 0, 0) \, \mathrm{Gauss}$. In this case, we need to apply Euler's method twice within each iteration of the integration. First, we calculate the new position of the electron given its instantaneous velocity:

$$\vec{x}_{i+1} = \vec{x}_{i} + h \vec{v}_{i}$$

and then we calculate the new velocity of the electron given its instantaneous acceleration:

$$\vec{v}_{i+1} = \vec{v}_{i} + h \vec{a}_{i} = \vec{v}_{i} + h \frac{e}{m_e} \left( \frac{\vec{v}_{i}}{c} \times \vec{B} \right)$$.

In our code, we will use a step size of 1e-6 seconds and integrate for a total of 0.1 second.

\subsection*{ The Runge-Kutta Method }

\section{ Monte Carlo Methods }

\subsection*{ Monte Carlo Sampling }

\subsection*{ Markov Chain Monte Carlo }

\end{document} <\latex>