""" They are separated by hashtag lines. Note that these require the odespy package to be installed. Open a Canopy terminal from the Tools menu in Canopy, and run pip install -e git+https://github.com/hplgit/odespy.git#egg=odespy """ #Solver Method Example import numpy import matplotlib.pyplot as plt import odespy from math import sqrt pi = 3.141592653589793 class Problem: def __init__(self, c, Theta): self.c, self.Theta = float(c), float(Theta) self.freq = sqrt(c) self.period = 2*pi/self.freq def f(self, u, t): theta, omega = u; c = self.c return [omega, -c*theta] problem = Problem(c=1, Theta=pi/4) solvers = [ odespy.ThetaRule(problem.f, theta=0), # Forward Euler odespy.ThetaRule(problem.f, theta=0.5), # Midpoint method odespy.ThetaRule(problem.f, theta=1), # Backward Euler odespy.RK4(problem.f), # Runge Kutta odespy.MidpointIter(problem.f, max_iter=2, eps_iter=0.01), odespy.LeapfrogFiltered(problem.f) ] Ns = [20,50,100,500] for q in range(1,5): N_per_period = Ns[q-1] T = 3*problem.period # final time legends = [] for solver in solvers: solver_name = str(solver) # short description of solver print solver_name solver.set_initial_condition([problem.Theta, 0]) N = N_per_period*problem.period time_points = numpy.linspace(0, T, N+1) u, t = solver.solve(time_points) theta = u[:,0] legends.append(solver_name) plt.subplot(2,2,q) plt.plot(t, theta) plt.hold('on') if q==2: plt.legend(legends) plt.show()