These examples are adapted from Sage wiki: https://wiki.sagemath.org/SageWiki
x,y = var('x,y')
The following two cells implement three integration methods: midpoint, trapezoidal, Simpson's rules. In all the three methods the interval $[-1,1]$ is partitioned into $n$ subintervals. The corresponding sums consist of n (midpoint) and n+1 terms (trapezoidal, Simpson's rule), respectively.
import scipy
import numpy
from scipy.special.orthogonal import p_roots, t_roots, u_roots
from scipy.integrate import quad, trapz, simps
from sage.ext.fast_eval import fast_float
from numpy import linspace
methods = {
'Midpoint': {'w': 1, 'xmin': -1, 'xmax': 1,
'func': lambda n: (linspace(-1r+1.0r/n,1r-1.0r/n,n),
numpy.array([1.0r]*n)*2.0r/n
),
'nterms': lambda n: n},
'Trapezoid': {'w': 1, 'xmin': -1, 'xmax': 1,
'func': lambda n: (linspace(-1r,1r,n+1),
numpy.array([1.0r]+[2.0r]*(n-1)+[1.0r])*1.0r/n
),
'nterms': lambda n: n+1},
'Simpson': {'w': 1, 'xmin': -1, 'xmax': 1,
'func': lambda n: (linspace(-1r,1r,n+1),
numpy.array([1.0r]+[4.0r,2.0r]*int((n-1.0r)/2.0r)+[4.0r,1.0r])*2.0r/(3.0r*n)
),
'nterms': lambda n: n}}
var("x")
def box(center, height, width2,**kwds):
return polygon([(center-width2,0),
(center+width2,0),(center+width2,height),(center-width2,height)],**kwds)
@interact
def weights(n=slider(2,30,2,default=10),
f=input_box(default=cos(x),type=SR),
show_method=["Trapezoid","Midpoint","Simpson"]):
ff = fast_float(f,'x')
method = methods[show_method]
xcoords,w = (method['func'])(int(n))
xmin = method['xmin']
xmax = method['xmax']
plot_min = max(xmin, -10)
plot_max = min(xmax, 10)
coords = zip(xcoords,w)
max_weight = max(w)
coords_scaled = zip(xcoords,w/max_weight)
f_graph = plot(f,plot_min,plot_max)
# boxes = sum(box(x,ff(x),(xmax-xmin)/n/2,rgbcolor=(0.5,0.5,0.5),alpha=0.3) for x,w in coords)
stems = sum(line([(x,0),(x,ff(x))],rgbcolor=(1-y,1-y,1-y),
thickness=2,markersize=6,alpha=y) for x,y in coords_scaled)
points = sum([point([(x,0),
(x,ff(x))],rgbcolor='black',pointsize=30) for x,_ in coords])
graph = stems+points+f_graph#+boxes
show(graph,xmin=plot_min,xmax=plot_max,aspect_ratio="auto")
approximation = sum([w*ff(x) for x,w in zip(xcoords,w)])
integral,integral_error = scipy.integrate.quad(ff, xmin, xmax)
x_val = linspace(-1, 1,n)
y_val = [*map(ff,x_val)]
trapezoid = integral-trapz(y_val, x_val)
simpson = integral-simps(y_val, x_val)
pretty_print(html(r"$$\sum_{i=1}^{i=%s} w_i f(x_i^*) \Delta x = %s\approx \int_{-1}^{1}%s \,dx = %s $$"%
( (method['nterms'])(int(n)), approximation, latex(f), integral)))
error_data = [integral-approximation, trapezoid, simpson, integral_error]
print("Errors for integration methods:\n")
print("Our implementation:\n%s\t %s\n\nBuilt-in:\nquad\t\t %s\ntrapz\t\t %s\nsimps\t\t %s" % (
show_method,integral-approximation,integral_error,trapezoid,simpson))
In the formula above, $w_i$ are the weights on values of $f$, which show up in the integration rules. For example, we have seen that the n-th sum in the midpoint rule is $$ M_n = \Delta x\left(f(\bar x_1) + f(\bar x_2)+\ldots+ f(\bar x_n) \right). $$ We interpret this expression as the sum $$ \sum_{i=1}^n w_i f(\bar x_i)\,\Delta x, $$ where all the $w_i$ are equal to 1. On the other hand, for the trapezoidal rule, the n-th sum is $$ M_n = \frac{\Delta x}{2} \left(f(x_0) + 2f(\bar x_1)+\ldots+ 2f(x_{n-1}) + f(x_{n}) \right), $$ so in this sum $$ \sum_{i=0}^n w_i f(x_i)\,\Delta x, $$ $w_i$ are equal to $$ \frac 12, \ 1, \ 1, \ 1, \ldots, 1,\ 1, \ \frac12. $$ (Just multiply the 1/2 in front with the coefficients at $f(x_i)$.)
The following cell implements Euler's method and compares the approximate solution obtained with it (blue) to the exact solution (red).
@interact
def euler_method(
y_exact_in = input_box('exp(x)', type = str, label = 'Exact solution = '),
y_prime_in = input_box('exp(x)', type = str, label = "y' = "),
start = input_box(0.0, label = 'x starting value: '),
startval = input_box(1.0, label = 'y starting value: '),
stop = input_box(3.0, label = 'x stopping value: '),
nsteps = slider([2^m for m in range(0,10)], default = 10, label = 'Number of steps: '),
):
y_exact = lambda x: eval(y_exact_in)
# dy/dx = F(x,y)
y_prime = lambda x,y: eval(y_prime_in)
stepsize = float((stop-start)/nsteps)
sol = [startval]+[0.0]*nsteps
xvals = [start]+[0.0]*nsteps
for step in range(nsteps):
sol[step+1] = sol[step] + stepsize*y_prime(xvals[step],sol[step])
xvals[step+1] = xvals[step] + stepsize
sol_max = max(sol + [find_local_maximum(y_exact,start,stop)[0]])
sol_min = min(sol + [find_local_minimum(y_exact,start,stop)[0]])
show(plot(y_exact(x),start,stop,rgbcolor=(1,0,0))+line([[xvals[index],sol[index]] for index in range(len(sol))]),xmin=start,xmax = stop, ymax = sol_max, ymin = sol_min)
In the next two cells, we plot the vector field for the vector with coordinates $(f,g)$, where both $f$ and $g$ depend on $(x,y)$. The vector field puts at every point $(x,y)$ of the plain vector $$ (f(x,y), g(x,y).$$ In the setting from the lecture, where we had $$ \frac{dy}{dx} = F(x,y), $$ this means that $$ F(x,y) = \frac{g(x,y)}{f(x,y)}. $$ That is, the slope of a vector is given by the ratio of second coordinate to the first one.
The blue curve connects the points $(x_k, y_k)$ for $k\geq 0$, obtained by the iterations of the Euler's method. By making the step size ($\Delta x$ in the lecture) sufficiently small, we get a good approximation to the solution of the equation, corresponding to the vector field.
from sage.ext.fast_eval import fast_float
@interact
def _(f = input_box(default=y), g=input_box(default=-x*y+x^3-x),
xmin=input_box(default=-1), xmax=input_box(default=1),
ymin=input_box(default=-1), ymax=input_box(default=1),
start_x=input_box(default=0.5), start_y=input_box(default=0.5),
step_size= slider([2^m for m in range(-10,1)], default =0.125, label = r'$\Delta x$ '),
steps=slider([2^m for m in range(0,11)], default =128, label = 'Number of steps: ') ):
ff = fast_float(f, 'x', 'y')
gg = fast_float(g, 'x', 'y')
steps = int(steps)
points = [ (start_x, start_y) ]
for i in range(steps):
xx, yy = points[-1]
try:
points.append( (xx + step_size * ff(xx,yy), yy + step_size * gg(xx,yy)) )
except (ValueError, ArithmeticError, TypeError):
break
starting_point = point(points[0], pointsize=50)
solution = line(points)
vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )
result = vector_field + starting_point + solution
pretty_print(
html(
r"$\displaystyle\frac{dy}{dx} = \frac gf = \frac{%s}{%s}$ <br> " % (latex(g),latex(y))
)
)
result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)