8 Week 8 Participation

Name: [Type Your Name Here]

To begin all assignments (whether participation or homework), please save a copy of this notebook to your Google Drive by clicking File -> Save a copy in Drive

import numpy as np
import plotly.graph_objects as go
from scipy import optimize

from ipywidgets import interactive
import matplotlib.pyplot as plt

Practice (By Hand)

For the following differential equation, compute the first 3 terms using Euler’s method and the Midpoint Method

\[ x' = - \frac{1}{2} x^2, \; x(0) = 6 \]

Euler’s method $x(1) = , x(2) = , x(3) = $

Midpoint method $x(1) = , x(2) = , x(3) = $

\[ x(0.5) = x(0) + 0.5 * x'(0) \\ = 6 + 0.5(-1/2 \; x(0)^2 ) \\ = 6 + 0.5 (-18)\\ = -3\\ x(1) = x(0) + 1 * x'(0.5) \\ = 6 + 1 (-1/2 * x(0.5)^2)\\ = 6 + 1(-9/2) \\ = 1.5 \]

We will begin by adding a docstring for the functions we have already built.

def euler(f,x0,t0,tmax,dt):
  N = int(np.floor((tmax-t0)/dt)+1)
  t = np.linspace(t0,tmax,N)
  x = np.zeros(len(t))
  x[0] = x0
  for n in range(len(t)-1):
    x[n+1] = x[n] + dt*f(t[n],x[n])
  return t, x

def midpoint1d(f,x0,t0,tmax,dt):
    N = int(np.floor((tmax-t0)/dt)+1)
    t = np.linspace(t0,tmax,N)
    x = np.zeros(len(t)) # build an array for the x values
    x[0] = x0 # build the initial condition

    for n in range(len(t)-1):
      xtemp = x[n] + dt/2 * f(t[n], x[n])
      x[n+1] = x[n] + dt*f(t[n] + dt/2, xtemp)
    return t, x

8.1 Beyond Euler’s and Midpoint Method

As a recap, we have attempted to solve ODEs using the following approach: * make a discrete approximation to the derivative and * step forward through time as a difference equation.

Other higher-order methods include Runge-Kutta 4 methods. Even though we will not cover this in the course, these methods could be considered for the next midterm.

8.2 Backward Euler’s Method

Returning to

\[ x' = - \frac{1}{2} x^2, \; x(0) = 6 \]

Calculate the $x(1)$ and $x(2)$ using Backward Euler’s Method when $h=1$.

$x(1) = ???, x(2) = ???$

8.2.1 Side-Bar - `scipy.optimize.fsolve()`

The scipy.optimize.fsolve() function finds the roots of a function, and returns the roots of the (non-linear) equations defined by func(x) = 0 given a starting estimate.

Example

\[ y^2 + y = 4 \]

We can plot this and see where the zeros are: https://www.desmos.com/calculator

Now, let’s use optimize.fsolve

f1 = lambda y: y**2+y-4
z = optimize.fsolve(f1,1)
print(z)
# start at a different starting point
z = optimize.fsolve(f1,-2)
print(z)

[1.56155281]
[-2.56155281]

Use this approach to solve the following problems:

$x^3 + 5x - 3 = \sin(x)$
$\sin^2 (x) - \cos(x) = -e^x$
$\frac{1}{2}\left(\left|x\right|+x\right)$ = 0

Now, we will use this tool to help with our backwardEuler1d() function.

Practice

Based on our discussion, complete the code below:

def backwardEuler1d(f, x0, t0, tmax, dt):
  t = np.arange(t0, tmax+dt, dt)
  x = np.zeros(len(t))
  x[0] = x0
  for n in range(len(t)-1):
    G = lambda X: X - x[n] - dt*f(t[n+1], X) # define this function
    # give the correct starting point for the solver below
    x[n+1] = optimize.fsolve(G, x[n])[0]

  return t,x

Practice

Test your code on the $x' = -1/2 x^2$ problem to verify your work.

fp = lambda t, x: -0.5 * x**2

8.3 Fitting ODE Models to Data

Newton’s Law of Cooling states that

\[ \frac{dT}{dt} = -k(T - T_{ambient}), \]

where $T$ is the temperature of some object (like a cup of coffee), $T_{ambient}$ is the temperature of the ambient environment, and $k$ is the proportionality constant that governs the rate of cooling. This is a classic differential equation with a well known solution. In the present situation we don’t want the analytic solution, but instead we will work with a numerical solution since we are thinking ahead to where the differential equation may be very hard to solve in future problems.

We also don’t want to just look at the data and guess an algebraic form for the function that best fits the data. That would be a trap! (why?)

The following data table gives the temperature (degrees F) at several times while a cup of tea cools on a table. The ambient temperature of the room is 65F.

import pandas as pd

URL1 = 'https://raw.githubusercontent.com/NumericalMethodsSullivan'
URL2 = '/NumericalMethodsSullivan.github.io/master/data/'
URL = URL1+URL2
data = pd.read_csv(URL+'Exercise5_newtoncooling.csv')
# Exercise5_newtoncooling.csv
#
data
data_np = np.array(data)
# or you can load the data directly with
# data = np.array([[0,160],[60,155],[180,145],[210,142],[600,120]])
fig = go.Figure()
fig.add_trace(go.Scatter(x=data['Time (sec)'], y=data['Temperature'], mode='markers'))

fig.update_layout(title="Tea Cooling", xaxis_title="Time (sec)", yaxis_title="Temperature")

Tambient = 65 # temp of room

# next we define our specific differential equation
f = lambda t, x, k: -k*(x - Tambient)
x0 = 160 # initial condition pulled from the data

def numericalSolution(k, dt):
    t0 = 0 # time where the data starts
    tmax = 660 # just beyond where the data ends
    t = np.arange(t0,tmax+dt,dt) # set up the times
    x = np.zeros(len(t))
    x[0] = x0
    for n in range(len(t)-1):
        # put the code necessary to build a good
        # numerical solver here be sure to account
        # for the parameter k in each of your function calls.
        # i.e., euler/midpoint
        x[n+1] = x[n] + dt*f(t[n], x[n], k)
    return t, x

data

	Time (sec)	Temperature
0	0	160
1	60	155
2	180	145
3	210	142
4	600	120

t = np.arange(0, 660 + 1, 1)

# this code block determines where the indices for the data match
indices = []
for j in range(len(t)):
    for k in range(len(data_np)):
        if t[j] == data_np[k,0]:
            indices.append(j)

indices

[0, 6, 18, 21, 60]

t[21]

def dataMatcher(k):
    # Next choose an appropriate value of dt.
    # Choosing dt so that values of time in the data fall within
    # the times for the numerical solution is typically a good
    # practice (but is not always possible).
    dt = ...
    # approximate solution
    t, x = numericalSolution(k, dt)
    # keep track of error at the data supplied times
    err = []
    counter = 0
    for n in indices:
       # calculate the error and save it
        err.append( (data_np[counter,1] - x[int(n)])**2 )
        counter += 1
    #print("For k=",k[0],",  SSRes=",np.sum(err)) # optional
    return np.sum(err)

import scipy.optimize as sp
# Choose an initial value of k and put it into the following code
# in place of the "???".  Note that we are sending a few parameters
# to the optimization tool.  Be sure to understand these options
# and take care that these options problem dependent and you will
# need to choose these again for the next new problem.
k_initial = ???
K = sp.minimize(dataMatcher, k_initial, options = {'maxiter': 5}, tol=1e-2)
print(K)
dt = 10
t1, x = numericalSolution(K.x[0], dt)
# plots
fig = go.Figure()
fig.add_trace(go.Scatter(x=data['Time (sec)'], y=data['Temperature'], name='Data', mode='markers'))
fig.add_trace(go.Scatter(x=t1, y=x, name='approximate solution'))
fig.update_layout(title="Tea Cooling", xaxis_title="Time (sec)", yaxis_title="Temperature")

  message: Maximum number of iterations has been exceeded.
  success: False
   status: 1
      fun: 1.8524481602331837
        x: [ 9.313e-04]
      nit: 5
      jac: [ 1.988e+02]
 hess_inv: [[ 3.191e-10]]
     nfev: 30
     njev: 10

RuntimeWarning:

overflow encountered in double_scalars

/usr/local/lib/python3.10/dist-packages/scipy/optimize/_numdiff.py:576: RuntimeWarning:

invalid value encountered in subtract

<ipython-input-28-6a25d2d05ddf>:14: RuntimeWarning:

overflow encountered in double_scalars

8.1 Beyond Euler’s and Midpoint Method

8.2 Backward Euler’s Method

8.2.1 Side-Bar - scipy.optimize.fsolve()

8.3 Fitting ODE Models to Data

8.2.1 Side-Bar - `scipy.optimize.fsolve()`