Question

1. Write a python code that uses the Runge Kutta Method method to approximate the solutions to each of the following initial-value problems and compare/plot the results to the actual values.

a) y′=te^(3t) − 2y, 0 < t < 1, y(0) = 0

with h = 0.5; actual solution y(t)=1/5te^(3t) − 1/25e^(3t) + 1/25e^(−2t).

- Use the Runge Kutta method to approximate/plot the solutions to each of the following initial-value

b) ?′=1+(?−?)2,2<?<3,?(2)=1y′=1+(t−y)2,2

c) ?′=1+??,1<?<1,?(1)=2y′=1+yt,1

Answer 1

# Some basic functional programming libraries
from functools import reduce
from operator import mul

def runge_kutta(f, t0, y0, h):
'''
Classical Runge-Kutta method for dy/dt = f(t, y), y(t0) = y0,
with step h, and the specified tolerance and max_steps.
This function is a generator which can give infinitely many points
of the estimated solution, in pairs (t[n], y[n]).
To get only finitely many values of
the solution we can for example do,
>>> from itertools import islice
>>> list(islice(runge_kutta(f, t0, h), n))
[(t[0], y[0]), (t[1], y[1]), ..., (t[n], y[n])]
and infact, we could define another function to do this like,
>>> runge_kutta_N = lambda f, t0, y0, h, N: list(islice(
... runge_kutta(f, t0, y0, h), N))
It would also be easy to change this function to take an extra
parameter N and then return a list of the first N, (t_n, y_n),
directly (by replacing the while loop with for n in range(N)).
Note also that whilst the question asks for a solution, this
function only returns an approximation of the solution at
certain points. We can turn use this to generate a continuous
function passing through the points specified using either of
interpolation methods specified lower down the file.
'''
# y and t represent y[n] and t[n] respectively at each stage
y = y0
t = t0

# Whilst it would be more elegant to write this recursively,
# in Python this would be very inefficient, and lead to errors when
# many iterations are required, as the language does not perform
# tail call optimisations as would be the case in languages such
# as C, Lisp, or Haskell.
#
# Instead we use a simple infinite loop, which will yield more values
# of the function indefinitely.
while True:
# Generate the next values of the solution y
yield t, y

# Values for weighted average (compare with Wikipedia)
k1 = f(t, y)
k2 = f(t + h/2, y + (h/2)*k1)
k3 = f(t + h/2, y + (h/2)*k2)
k4 = f(t + h/2, y + h*k3)

# Calculate the new value of y and t as per the Runge-Kutta method
y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
t = t + h

def linear_interpolation(tys):
'''Takes a list of (t, y(t)) values (presumed to be in increasing
order of t), and interpolates to a piecewise linear function
passing through each point.
Heavily inspired by: http://en.wikipedia.org/wiki/Linear_interpolation
'''
# The t values
ts = [t for t, y in tys]
# The y values
ys = [y for t, y in tys]
# The number of values
n = len(tys) # = len(ts) = len(ys)

if n < 2:
raise ValueError("Not enough points to interpolate!")

def z(x):
# If x is beyond the lower end of the data, we extrapolate
# based on the first two values.
if x <= min(ts):
i = 0
return ys[0] + (ys[1] - ys[0]) * (x - ts[0])/(ts[1] - ts[0])
# If x is beyond the upper end of the data, we extrapolate
# based on the last two values.
elif x >= max(ts):
i = n-2
else:
# Find the interval i, that x lies within
i = next(i for i in range(n - 1) if ts[i] <= x <= ts[i+1])

# Interpolate within our chosen interval
return ys[i] + (ys[i+1] - ys[i]) * (x - ts[i])/(ts[i+1] - ts[i])

# Return a continous function z which interpolates between the given values
return z

# A concise definition of the product of a list/iterable of numbers
# (sum is built in but product is not)
product = lambda xs: reduce(mul, xs, 1)

def polynomial_interpolation(tys):
'''
Takes a list of n (t, y(t)) values (presumed to be in increasing
order of t), and interpolates to a nth order polynomial
passing through each point.
Uses the Lagrange polynomial:
- http://en.wikipedia.org/wiki/Lagrange_polynomial
- http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html
Note that this implementation could easily be made much faster
at the expense of clarity (and similarity to the Wikipedia
article).
'''
# The t values
ts = [t for t, y in tys]
# The y values
ys = [y for t, y in tys]
# The number of values
n = len(tys) # = len(ts) = len(ys)

# Return the interpolation function
return lambda x: sum(ys[j]*product((x - ts[k])/(ts[j] - ts[k])
for k in range(n) if j != k) for j in range(n))

if __name__ == '__main__':
# A demonstration of the method for,
# dy/dt = (1/3)*y + e^t*y^4, y(0) = 1
# (chosen to test the method, as an explict solution does exist)
print('dy/dt = (1/3)*y + e^t*y^4, y(0) = 1')
from itertools import islice
from math import exp
# Evaluate the solution at some points
tys = list(islice(runge_kutta(
lambda t, y: (1/3)*y + exp(t)*y**4, # f(t, y)
0, # t0
1, # y0
0.001,
), 200))
# Print some values
print('\n'.join('y({:.3f}) = {}'.format(t, y) for t, y in tys))

# The explicit solution of the equation
z = lambda t: ((-2)**(1/3)*exp(t/3))/(3*exp(2*t) - 5)**(1/3)
# Calculate maximum absolute error
e_abs_max = max(abs(y - z(t)) for t, y in tys)
# Calculate maximum relative error
e_rel_max = max(abs((y - z(t))/z(t)) for t, y in tys)
# Print the errors
print('Maximum absolute error: ', e_abs_max)
print('Maximum relative error:', e_rel_max)

# To check the interpolation is not completely messed up, check it
# agrees on the actual data points
z_linear = linear_interpolation(tys)
z_polynomial = polynomial_interpolation(tys)
print('Linear interpolation agrees?',
'Yes' if all(y == z_linear(t) for t, y in tys) else 'No')
print('Polynomial interpolation agrees?',
'Yes' if all(y == z_polynomial(t) for t, y in tys) else 'No')

try:
# Plot some graphs to compare the results, this time with only
# 20 points to compare interpolation methods more clearly
import matplotlib.pyplot as plt

# Evaluate the solution at some points
tys2 = list(islice(runge_kutta(
lambda t, y: (1/3)*y + exp(t)*y**4, # f(t, y)
0, # t0
1, # y0
0.02,
), 10))
ts = [t for t, y in tys2]
ys = [y for t, y in tys2]
z_linear_2 = linear_interpolation(tys2)
z_polynomial_2 = polynomial_interpolation(tys2)
ts_dense = [n*0.001 for n in range(200)]
plt.plot(ts_dense, [z_linear_2(t) for t in ts_dense], 'g',
label='Linear')
plt.plot(ts_dense, [z_polynomial_2(t) for t in ts_dense], 'r',
label='Polynomial')
plt.plot(ts_dense, [z(t) for t in ts_dense], 'y',
label='Exact solution')
plt.plot(ts, ys, 'b^', label='Runge-Kutta')
plt.legend()
plt.ylabel("y")
plt.xlabel("t")
plt.show()
except ImportError:
print("Matplotlib required for graphs")

Is also Try:

import numpy as np
from math import sqrt

import matplotlib.pyplot as plt

# Parameters
mr = 2
m = 1
c = 5
k = 36
F = 3

# Equations:

def V(u,t):
x1, x2, v1, v2 = u
dx = x2-x1;
dv = v2-v1;

return np.array([ v1, v2, (c*dv+k*dx)/mr, -(F+c*dv+k*dx)/m ])

def rk4(f, u0, t0, tf , n):
t = np.linspace(t0, tf, n+1)
u = np.array((n+1)*[u0])
h = t[1]-t[0]
for i in range(n):
k1 = h * f(u[i], t[i])
k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
k4 = h * f(u[i] + k3, t[i] + h)
u[i+1] = u[i] + (k1 + 2*(k2 + k3) + k4) / 6
return u, t

u, t = rk4(V, np.array([0., 0., 1., 1.]) , 0. , 1. , 10)
x1, x2, v1, v2 = u.T
plt.plot(t, x1, t, x2)
plt.grid('on')
plt.show()