Update documentation

MAT301_DifferentialEquationSolving1.html

@@ -262,6 +262,11 @@ <h1 class="site-logo" id="site-title">MAT301 Numerical Methods</h1>
    Euler’s method
   </a>
  </li>
+ <li class="toc-h2 nav-item toc-entry">
+  <a class="reference internal nav-link" href="#example-9-1">
+   Example 9.1
+  </a>
+ </li>
  <li class="toc-h2 nav-item toc-entry">
   <a class="reference internal nav-link" href="#midpoint-method">
    Midpoint method
@@ -316,6 +321,11 @@ <h2> Contents </h2>
    Euler’s method
   </a>
  </li>
+ <li class="toc-h2 nav-item toc-entry">
+  <a class="reference internal nav-link" href="#example-9-1">
+   Example 9.1
+  </a>
+ </li>
  <li class="toc-h2 nav-item toc-entry">
   <a class="reference internal nav-link" href="#midpoint-method">
    Midpoint method
@@ -344,7 +354,7 @@ <h2> Contents </h2>
 
   <div class="tex2jax_ignore mathjax_ignore section" id="differential-equation-solving-1">
 <h1>Differential Equation Solving 1<a class="headerlink" href="#differential-equation-solving-1" title="Permalink to this headline">¶</a></h1>
-<p>Notebook to implement Euler, Midpoint and Euler Trapezium methods to solve differential equations. Designed to support Abertay University undergrad course MAT301 (JT, 2022).</p>
+<p>Many physical laws governing natural behaviour take the form of differential equations. These are often complicated, coupled equations where the behaviour of quantities in space and time depend on others. While some have analytical solutions, often it is quicker and easier to apply numerical tools to obtain approximate solutions.</p>
 <div class="section" id="motivation">
 <h2>Motivation<a class="headerlink" href="#motivation" title="Permalink to this headline">¶</a></h2>
 <p>A <strong>differential equation</strong> is a relationship between a function, <span class="math notranslate nohighlight">\(f(x)\)</span>, its independent variable, <span class="math notranslate nohighlight">\(x\)</span>, and any number of its derivatives. An <strong>ordinary differential equation</strong> or <strong>ODE</strong> is a differential equation where the independent variable, and therefore also the derivatives, is in one dimension. For our purposes, we can assume that an ODE can be written</p>
@@ -395,8 +405,11 @@ <h2>Euler’s method<a class="headerlink" href="#euler-s-method" title="Permalin
 \[
 S(t_{j+1}) = S(t_j) + hF(t_j, S(t_j)).
 \]</div>
-<p>This formula is called the <strong>Explicit Euler Formula</strong>, and allows us to approximate the state at <span class="math notranslate nohighlight">\(S(t_{j+1})\)</span> given the state at <span class="math notranslate nohighlight">\(S(t_j)\)</span>. Starting from a given initial value of <span class="math notranslate nohighlight">\(S_0 = S(t_0)\)</span>, we can use this formula to integrate the states up to <span class="math notranslate nohighlight">\(S(t_f)\)</span>; these <span class="math notranslate nohighlight">\(S(t)\)</span> values are then an approximation for the solution of the differential equation. The Explicit Euler formula is the simplest and most intuitive method for solving initial value problems. At any state <span class="math notranslate nohighlight">\((t_j, S(t_j))\)</span> it uses <span class="math notranslate nohighlight">\(F\)</span> at that state to “point” toward the next state and then moves in that direction a distance of <span class="math notranslate nohighlight">\(h\)</span>. Although there are more sophisticated and accurate methods for solving these problems, they all have the same fundamental structure. As such, we enumerate explicitly the steps for solving an initial value problem using the Explicit Euler formula.</p>
-<img src="https://pythonnumericalmethods.berkeley.edu/_images/22.03.01-Euler-method-illustration.png" alt="explicit Euler" title="The illustration of the explicit Euler method." width="500"/><p>Lets go ahead and see this in action, using the equation and question posed in from Example 9.1 in the lectures:</p>
+<p>This formula is called the <strong>Explicit Euler Formula</strong>, and allows us to approximate the state at <span class="math notranslate nohighlight">\(S(t_{j+1})\)</span> given the state at <span class="math notranslate nohighlight">\(S(t_j)\)</span>. Starting from a given initial value of <span class="math notranslate nohighlight">\(S_0 = S(t_0)\)</span>, we can use this formula to integrate the states up to <span class="math notranslate nohighlight">\(S(t_f)\)</span>; these <span class="math notranslate nohighlight">\(S(t)\)</span> values are then an approximation for the solution of the differential equation. The Euler formula is the simplest and most intuitive method for solving initial value problems. At any state <span class="math notranslate nohighlight">\((t_j, S(t_j))\)</span> it uses <span class="math notranslate nohighlight">\(F\)</span> at that state to “point” toward the next state and then moves in that direction a distance of <span class="math notranslate nohighlight">\(h\)</span>. Although there are more sophisticated and accurate methods for solving these problems, they all have the same fundamental structure. As such, we enumerate explicitly the steps for solving an initial value problem using the Explicit Euler formula.</p>
+<img src="https://pythonnumericalmethods.berkeley.edu/_images/22.03.01-Euler-method-illustration.png" alt="explicit Euler" title="The illustration of the explicit Euler method." width="500"/></div>
+<div class="section" id="example-9-1">
+<h2>Example 9.1<a class="headerlink" href="#example-9-1" title="Permalink to this headline">¶</a></h2>
+<p>To see this in action, we will use the equation and question posed in from Example 9.1 in the lectures:</p>
 <div class="math notranslate nohighlight">
 \[
 \frac{{\rm{d}}y}{{\rm{d}}t}=xy,
@@ -488,7 +501,11 @@ <h2>Euler’s method<a class="headerlink" href="#euler-s-method" title="Permalin
 </div>
 </div>
 <div class="cell_output docutils container">
-<img alt="_images/MAT301_DifferentialEquationSolving1_17_0.png" src="_images/MAT301_DifferentialEquationSolving1_17_0.png" />
+<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/9y/379gbxjx3255qz9rg64kh7mm0000gn/T/ipykernel_52116/2040643597.py:1: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as &#39;seaborn-v0_8-&lt;style&gt;&#39;. Alternatively, directly use the seaborn API instead.
+  plt.style.use(&#39;seaborn-poster&#39;)
+</pre></div>
+</div>
+<img alt="_images/MAT301_DifferentialEquationSolving1_17_1.png" src="_images/MAT301_DifferentialEquationSolving1_17_1.png" />
 </div>
 </div>
 <p>So we see that our approximation displays similar behaviour, but doesn’t follow the trend exactly; in fact the graphs appear to be diverging as <span class="math notranslate nohighlight">\(x\)</span> increases. What if we make the step size ten times smaller?</p>
@@ -516,20 +533,25 @@ <h2>Euler’s method<a class="headerlink" href="#euler-s-method" title="Permalin
 </div>
 </div>
 <div class="cell_output docutils container">
-<img alt="_images/MAT301_DifferentialEquationSolving1_19_0.png" src="_images/MAT301_DifferentialEquationSolving1_19_0.png" />
+<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/9y/379gbxjx3255qz9rg64kh7mm0000gn/T/ipykernel_52116/4080072740.py:8: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as &#39;seaborn-v0_8-&lt;style&gt;&#39;. Alternatively, directly use the seaborn API instead.
+  plt.style.use(&#39;seaborn-poster&#39;)
+</pre></div>
+</div>
+<img alt="_images/MAT301_DifferentialEquationSolving1_19_1.png" src="_images/MAT301_DifferentialEquationSolving1_19_1.png" />
 </div>
 </div>
 <p>As expected, cranking up the number of steps improves the accuracy of the solution. However, in practice, we often have a finite amount of numerical resources (e.g. CPU speed, number of nodes, memory size), so we can’t always rely on the ability to simply decrease <span class="math notranslate nohighlight">\(h\)</span>. In those cases, we have to be smarter, and apply other techniques.</p>
 </div>
 <div class="section" id="midpoint-method">
 <h2>Midpoint method<a class="headerlink" href="#midpoint-method" title="Permalink to this headline">¶</a></h2>
-<p>In the lectures, we were introduced to a more accurate scheme, the Midpoint method. This essentially estimates the gradient of the tangent at the midpoint between steps <span class="math notranslate nohighlight">\(x_n+\frac{1}{2}h\)</span>, which lies halfway between <span class="math notranslate nohighlight">\(x_n\)</span> and <span class="math notranslate nohighlight">\(x_{n+1}\)</span>. Using Euler’s method, we then determine \begin{equation}y_{n+1}=y_n+hk_2,\end{equation}where
-$<span class="math notranslate nohighlight">\(
+<p>In the lectures, we were introduced to a more accurate scheme, the Midpoint method. This essentially estimates the gradient of the tangent at the midpoint between steps <span class="math notranslate nohighlight">\(x_n+\frac{1}{2}h\)</span>, which lies halfway between <span class="math notranslate nohighlight">\(x_n\)</span> and <span class="math notranslate nohighlight">\(x_{n+1}\)</span>. Using Euler’s method, we then determine <span class="math notranslate nohighlight">\(y_{n+1}=y_n+hk_2,\)</span> where</p>
+<div class="math notranslate nohighlight">
+\[\begin{split}
 \begin{align}
-k_1&amp;=f(x_n,y_n) \\
+k_1&amp;=f(x_n,y_n), \\
 k_2&amp;=f(x_n+\frac{1}{2}h, y_n+\frac{1}{2}hk_1).
 \end{align}
-\)</span>$</p>
+\end{split}\]</div>
 <p>This looks a bit weird, but all this means is that we start to evaluate the function not only at the grid point values, but also now at different locations in <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span>.</p>
 <p>Let’s go ahead and implement this, in exactly the same manner as the previous method:</p>
 <div class="cell docutils container">
@@ -555,7 +577,11 @@ <h2>Midpoint method<a class="headerlink" href="#midpoint-method" title="Permalin
 </div>
 </div>
 <div class="cell_output docutils container">
-<img alt="_images/MAT301_DifferentialEquationSolving1_23_0.png" src="_images/MAT301_DifferentialEquationSolving1_23_0.png" />
+<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/9y/379gbxjx3255qz9rg64kh7mm0000gn/T/ipykernel_52116/2590250041.py:7: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as &#39;seaborn-v0_8-&lt;style&gt;&#39;. Alternatively, directly use the seaborn API instead.
+  plt.style.use(&#39;seaborn-poster&#39;)
+</pre></div>
+</div>
+<img alt="_images/MAT301_DifferentialEquationSolving1_23_1.png" src="_images/MAT301_DifferentialEquationSolving1_23_1.png" />
 </div>
 </div>
 <p>We quickly see that the midpoint method is more accurate than the Euler method for the same step size. The trade-off is the additional calculation required to evaluate <span class="math notranslate nohighlight">\(k_1\)</span> and <span class="math notranslate nohighlight">\(k_2\)</span>.</p>
@@ -571,12 +597,13 @@ <h3>Modified Euler (Euler-Trapezium Method)<a class="headerlink" href="#modified
 <p>where <span class="math notranslate nohighlight">\(y_{n+1}^p\)</span> is an initial prediction. We then use the derivative prediction, <span class="math notranslate nohighlight">\(y_{n+1}'^p\)</span>, from</p>
 <div class="math notranslate nohighlight">
 \[
-y_{n+1}'^p=f\left(x_{n+1}, y_{n+1}^p\right)
+y_{n+1}'^p=f\left(x_{n+1}, y_{n+1}^p\right),
 \]</div>
-<p>to correct the value
-$<span class="math notranslate nohighlight">\(
+<p>to correct the value</p>
+<div class="math notranslate nohighlight">
+\[
 y_{n+1}^c=y_n+\frac{h}{2}\left(y_n'+y_{n+1}'^p\right).
-\)</span>$</p>
+\]</div>
 <p>Once again, we can apply this solution to the ODE problem from earlier in the following way:</p>
 <div class="cell docutils container">
 <div class="cell_input docutils container">
@@ -601,7 +628,11 @@ <h3>Modified Euler (Euler-Trapezium Method)<a class="headerlink" href="#modified
 </div>
 </div>
 <div class="cell_output docutils container">
-<img alt="_images/MAT301_DifferentialEquationSolving1_27_0.png" src="_images/MAT301_DifferentialEquationSolving1_27_0.png" />
+<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/var/folders/9y/379gbxjx3255qz9rg64kh7mm0000gn/T/ipykernel_52116/1612707178.py:7: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as &#39;seaborn-v0_8-&lt;style&gt;&#39;. Alternatively, directly use the seaborn API instead.
+  plt.style.use(&#39;seaborn-poster&#39;)
+</pre></div>
+</div>
+<img alt="_images/MAT301_DifferentialEquationSolving1_27_1.png" src="_images/MAT301_DifferentialEquationSolving1_27_1.png" />
 </div>
 </div>
 <p>Once again, the Modified Euler does a much better job than the original Euler method, and (for this ODE) is comparable to the Midpoint method in accuracy.</p>

MAT301_DifferentialEquationSolving2.html

@@ -538,11 +538,12 @@ <h2>4th Order Runge Kutta Algorithm<a class="headerlink" href="#th-order-runge-k
 <p>To illustrate how this method can be implemented numerically, we will apply it to three initial value problems. These happen to be models of population growth over time, but could in reality be for any type of ordinary differential equation (ODE, or indeed for a system of several coupled ODEs). I have selected these equations because they are broadly related, but gradually increase in complexity. The equations which we will investigate are the following:</p>
 <div class="section" id="linear-population-equation">
 <h3>1. Linear Population Equation<a class="headerlink" href="#linear-population-equation" title="Permalink to this headline">¶</a></h3>
-<p>Our first ODE will be linear:
-$<span class="math notranslate nohighlight">\(
+<p>Our first ODE will be linear:</p>
+<div class="math notranslate nohighlight">
+\[
 y^{'}=0.1y, \ \ (2000 \leq t \leq 2020), 
-\)</span>$
-with the initial condition,</p>
+\]</div>
+<p>with the initial condition,</p>
 <div class="math notranslate nohighlight">
 \[
 y(2000)=6.
@@ -600,19 +601,21 @@ <h2>Discrete Interval<a class="headerlink" href="#discrete-interval" title="Perm
 \[
 h=\frac{b-a}{N}.
 \]</div>
-<p>Here, the interval is <span class="math notranslate nohighlight">\(2000\leq t \leq 2020,\)</span>
-$<span class="math notranslate nohighlight">\(
+<p>Here, the interval is <span class="math notranslate nohighlight">\(2000\leq t \leq 2020,\)</span></p>
+<div class="math notranslate nohighlight">
+\[
 h=\frac{2020-2000}{200}=0.1.
-\)</span>$</p>
+\]</div>
 <p>This gives the 201 discrete points:</p>
 <div class="math notranslate nohighlight">
 \[ 
 t_0=2000, \ t_1=2000.1, \ ... t_{200}=2020, 
 \]</div>
-<p>which we can generalise to
-$<span class="math notranslate nohighlight">\( 
+<p>which we can generalise to</p>
+<div class="math notranslate nohighlight">
+\[ 
 t_n=2000+0.1n, \ \ \ n=0,1,...,200.
-\)</span>$</p>
+\]</div>
 <p>Let’s go ahead and set this range of <span class="math notranslate nohighlight">\(t\)</span> values as a list in Python:</p>
 <div class="cell docutils container">
 <div class="cell_input docutils container">
@@ -668,14 +671,15 @@ <h3>Tailoring our 4th Order Runge Kutta scheme<a class="headerlink" href="#tailo
 \[
 k_1=f\left(x_n,y_n\right)=0.1 y_n.
 \]</div>
-<p>Subsequent coefficients are more complicated to evaluate mathematically, e.g. <span class="math notranslate nohighlight">\(k_2\)</span>:
-$<span class="math notranslate nohighlight">\(
+<p>Subsequent coefficients are more complicated to evaluate mathematically, e.g. <span class="math notranslate nohighlight">\(k_2\)</span>:</p>
+<div class="math notranslate nohighlight">
+\[\begin{split}
 \begin{align}
 k_2&amp;=f\left(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right)=0.1\left(y_n+0.05k_1\right) \\
 &amp;=0.1\left(y_n+0.05(0.1y_n)\right)=0.1y_n(1+0.005) \\
 &amp;=0.1y_n(1.005)=0.1005y_n.
 \end{align}
-\)</span>$</p>
+\end{split}\]</div>
 <p>This is a lot of mathematical effort to then plug into our calculators!</p>
 <p>The beauty of this method relies on making the computer to do the hard work for us! Now that we’ve defined a function to spit out values of <span class="math notranslate nohighlight">\(f(x,y)\)</span>, we can also evaluate <em>shifted values</em> of <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> as input into that function to evaluate <span class="math notranslate nohighlight">\(k_1-k_4\)</span>, and hence the next iteration of the solution, <span class="math notranslate nohighlight">\(y_{n+1}\)</span>:</p>
 <div class="cell docutils container">
@@ -908,14 +912,16 @@ <h3>Results and Graph<a class="headerlink" href="#results-and-graph" title="Perm
 \[
 \int{\frac{\rm{d}x}{x\left(a+bx\right)}}=\frac{1}{a}\ln{\left|\frac{a+bx}{x}\right|}+c,
 \]</div>
-<p>we can therefore find:
-$<span class="math notranslate nohighlight">\(
+<p>we can therefore find:</p>
+<div class="math notranslate nohighlight">
+\[
 y(t)=\left({5e^{\dfrac{-(t+c)}{5}}+0.05}\right)^{-1}, 
-\)</span>$</p>
-<p>with a contant of integration
-$<span class="math notranslate nohighlight">\(
+\]</div>
+<p>with a contant of integration</p>
+<div class="math notranslate nohighlight">
+\[
 c=-5\ln{\left(\frac{0.14\exp{(400)}}{6}\right)}.
-\)</span>$</p>
+\]</div>
 <p>The plot below again shows the RK4 numerical approximation, <span class="math notranslate nohighlight">\(y_n\)</span> (circles) and the analytical solution to the non-linear population equation:</p>
 <div class="cell docutils container">
 <div class="cell_input docutils container">
@@ -1058,16 +1064,18 @@ <h2>3. Non-Linear Population Equation with an oscillation<a class="headerlink" h
 \[ 
 y^{'}=0.2y-0.01y^2+\sin(2\pi t), \ \ (2000 \leq t \leq 2020), 
 \]</div>
-<p>with the initial condition,
-$<span class="math notranslate nohighlight">\(
+<p>with the initial condition,</p>
+<div class="math notranslate nohighlight">
+\[
 y(2000)=6.
-\)</span>$</p>
+\]</div>
 <div class="section" id="id3">
 <h3>Implementing RK4<a class="headerlink" href="#id3" title="Permalink to this headline">¶</a></h3>
-<p>As with the previous cases, lets express the right hand side of the ODE as a function
-$<span class="math notranslate nohighlight">\(
+<p>As with the previous cases, lets express the right hand side of the ODE as a function</p>
+<div class="math notranslate nohighlight">
+\[
 f_3(x,y)=0.2y-0.01y^2+\sin(2\pi x),
-\)</span>$</p>
+\]</div>
 <p>Once again, we could mathematically calculate expressions for <span class="math notranslate nohighlight">\(k_1\rightarrow k_4\)</span>.
 This time though, you can see that <span class="math notranslate nohighlight">\(f_3(x,y)\)</span> now contains terms that depend on <span class="math notranslate nohighlight">\(y\)</span> <strong>and</strong> <span class="math notranslate nohighlight">\(t\)</span>!</p>
 <p>Now that we’re more familiar with how to implement RK4 numerically, we’ll simply move on to define our function representing <span class="math notranslate nohighlight">\(f_3(x,y)\)</span> and let the code evaluate the coefficients:</p>