# Euler#

## 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.