Example 5: Newton's equation
Let's use Heun's method to solve Newton's equation from the section Introduction, y′(x)=1−3x+y+x2+xy, y(0)=0 with analytical solution y(x)=3√2πe⋅exp(x(1+x2))⋅[erf(√22(1+x))−erf(√22)]+4⋅[1−exp(x(1+x2))]−x
Here we have f(x,y)=1−3x+y+x2+xy=1+x(x−3)+(1+x)y
The following program NewtonHeun.py solves this problem using Heun's method, and the resulting figure is shown in Figure 7.
# chapter1/programs_and_modules/NewtonHeun.py
# Program Newton
# Computes the solution of Newton's 1st order equation (1671):
# dy/dx = 1-3*x + y + x^2 +x*y , y(0) = 0
# using Heun's method.
import numpy as np
xend = 2
dx = 0.1
steps = np.int(np.round(xend/dx, 0)) + 1
y, x = np.zeros((steps,1), float), np.zeros((steps,1), float)
y[0], x[0] = 0.0, 0.0
for n in range(0,steps-1):
x[n+1] = (n+1)*dx
xn = x[n]
fn = 1 + xn*(xn-3) + y[n]*(1+xn)
yp = y[n] + dx*fn
xnp1 = x[n+1]
fnp1 = 1 + xnp1*(xnp1-3) + yp*(1+xnp1)
y[n+1] = y[n] + 0.5*dx*(fn+fnp1)
# Analytical solution
from scipy.special import erf
a = np.sqrt(2)/2
t1 = np.exp(x*(1+ x/2))
t2 = erf((1+x)*a)-erf(a)
ya = 3*np.sqrt(2*np.pi*np.exp(1))*t1*t2 + 4*(1-t1)-x
# plotting
import matplotlib.pylab as py
py.plot(x, y, '-b.', x, ya, '-g.')
py.xlabel('x')
py.ylabel('y')
font = {'size' : 16}
py.rc('font', **font)
py.title('Solution to Newton\'s equation')
py.legend(['Heun', 'Analytical'], loc='best', frameon=False)
py.grid()
py.savefig('newton_heun.png', transparent=True)
py.show()
Figure 7: Velocity of falling sphere using Euler's and Heun's methods.
