Euler#

panelholoviews
Published: November 16, 2018 · Updated: October 31, 2023


Panel dashboard illustrating Euler’s Method#

Euler’s method is a numerical integration technique for solving ordinary differential equations. Specifically, given \(y' = f( x, y )\) and a starting point \((x_0, y_0)\), Euler’s method provides approximate values of \(y\) for \(x>x_0\). To explore how this method works, this notebook from Egbert Ammicht sets up a dashboard using Panel, HoloViews, and Bokeh. You will need to have a live Python process running (not just an Anaconda.org notebook or web page), which you can set up using conda install -c pyviz panel holoviews bokeh.

import param, panel as pn, holoviews as hv, numpy as np, pandas as pd
hv.extension('bokeh')
pn.extension()
#print(pn.__name__, pn.__version__, hv.__name__, hv.__version__ )

Vector field background#

To start, let’s consider the linear approximation \(y' \approx \frac{y(x+h)-y(x)}{h}\), where \(h\) is a small value controlling the integration step size. We can evaluate this equation on a set of \(x\) grid points \(x_n = x_0 + n h\), to obtain a series of \(y\) values using \(y_{n+1} = y_n + h f( x_n, y_n)\). We can use this approach to draw a vector field showing the direction of the integral at each location in the \(x,y\) space, overlaid on a grayscale plot of the same data:

def background(func, size=(500,500)):
    """
    Given the ODE y'=f(x,y),
    
       bg,vec,xaxis_line,yaxis_line = background()
    
    returns a grayscale image of the slopes y',
            a vector field representation of the slopes, and
            a set of axis lines for -5<x<5, -5<y<5
    """
    
    # compute the data
    vals = np.linspace( -5, 5, num=150 )
    X, Y = np.meshgrid(vals, vals)

    clines  = func(X,Y)                    # f(x,y)
    theta   = np.arctan2( clines, 1)       # angle of the slope y' at x,y

    # Obtain the vector field (subsample the grid)
    h,w=size
    vf_opts = dict(size_index=3, height=h, width=w, xticks=9, yticks=9, alpha=0.3, muted_alpha=0.05)
    vec_field = hv.VectorField((vals[::3],vals[::3], theta[::3,::3],0*clines[::3,::3]+1), 
                               label='vector_field' ).options(**vf_opts)
    
    # Normalize the given array so that it can be used with the RGB element's alpha channel   
    def norm(arr):
        arr = (arr-arr.min())
        return arr/arr.max()

    normXY    = norm(clines)
    img_field = hv.RGB( (vals, vals, normXY, normXY, normXY, 0*clines+0.1), vdims=['R','G','B','A'] )\
                .options(width=size[0], height=size[1], shared_axes=False)

    # finally, we add the axes as VLine, HLine and return an array of the plot Elements
    hv_opts = dict( color='k', alpha=0.8, line_width=1.5)
    return [img_field,vec_field, hv.HLine(0).options(**hv_opts),
                hv.VLine(0).options(**hv_opts)]

# Test it:
hv.Overlay(background(lambda x,y:   np.sin(x*y), size=(400,400))).options(show_legend=False).relabel("y' = sin(x y)"   ) +\
hv.Overlay(background(lambda x,y: x*np.sin(5*y), size=(400,400))).options(show_legend=False).relabel("y' = x sin(5 y)" )
WARNING:param.VectorFieldPlot: The `size_index` parameter is deprecated in favor of size style mapping, e.g. `size=dim('size')**2`.
WARNING:param.VectorFieldPlot: The `size_index` parameter is deprecated in favor of size style mapping, e.g. `size=dim('size')**2`.

Euler Integration Curves#

We can now overlay a series of lines approximating the integral for various step sizes at a given starting point, to show how the step size affects the accuracy.

def euler_step(x,y,h,func):
    """x <- x +h, y_<- y + h f(x,y)"""
    hs = h * func(x,y)
    x = x + h
    y = y + hs
    return x,y

def euler_table(x0,y0,n,h,func):
    """
    Compute up to n euler steps with step size h for  y' = f(x,y) starting from (x0,y0),
    returning the results in an hv.Table
    """
    xl = [x0]; yl=[y0]
    for i in range(n):
        x0,y0 = euler_step(x0,y0,h,func)
        xl.append(x0);yl.append(y0)
        if np.abs(x0) > 5. or np.abs(y0) > 5. : break   # we ran off the grid
    return hv.Table(pd.DataFrame(dict(x=xl,y=yl)), kdims=['x'],vdims=['y'])

line_colors = hv.Cycle(["Red","Orange","LightGreen","Green"])

def euler_curve(x0,y0,n,h,func):
    """
    Compute up to n euler steps with step size h for  y' = f(x,y) starting from (x0,y0)
    return the results in an hv.Curve
    """
    return euler_table(x0,y0,n,h,func).to.curve( label= 'h=%6.3f'%h).options(color=line_colors)

def append_euler_plots( l, start, func, n=10000, h=[.5,.2,.01,.0011] ):
    for hi in h: l.append( euler_curve(*start, n, hi, func) )
    return l
# Example functions

funcs = {
    "y' = sin(xy)"           : lambda x,y: np.sin(x*y),
    "y' = sin(x) sin(y)"      : lambda x,y: np.sin(x)*np.sin(y),
    "y' = cos(x)"             : lambda x,y: np.cos(x),
    "y' = exp(x /(x**2 + 1))" : lambda x,y: np.exp(-x/( x**2 + 1)),
    "y' = x**2 sin(5y)"      : lambda x,y: x**2*np.sin(5*y),
    "y' = x sin(5y)"         : lambda x,y: x   *np.sin(5*y),
    "y' = x tan(y)"           : lambda x,y: x   *np.tan(y),
    "y' = x / cosh(y)"        : lambda x,y: x/np.cosh(y),
}

f1_key   = "y' = sin(xy)"
l1       = background(funcs[f1_key])
append_euler_plots(l1, (-5,np.pi/4.75),funcs[f1_key] )

f2_key   = "y' = sin(x) sin(y)"
l2       = background(funcs[f2_key])
append_euler_plots(l2, (-5, np.pi/4.75),funcs[f2_key] )
append_euler_plots(l2, (-5,-np.pi/4.75),funcs[f2_key] )

# We need to call redim in case some curve overshot the grid  (might instead use apply_ranges=False to Curves)
pos_opts = dict(legend_position='right', toolbar='above',width=450,height=350)
hv.Overlay(l1).redim.range(x=(-5,5),y=(-5,5)).options(**pos_opts).relabel(f1_key) +\
hv.Overlay(l2).redim.range(x=(-5,5),y=(-5,5)).options(**pos_opts).relabel(f2_key)
WARNING:param.VectorFieldPlot: The `size_index` parameter is deprecated in favor of size style mapping, e.g. `size=dim('size')**2`.
WARNING:param.VectorFieldPlot: The `size_index` parameter is deprecated in favor of size style mapping, e.g. `size=dim('size')**2`.