I am trying to solve a Differential Equation with 4th Order Runge - Kutta method in Python3. The function is deliberately made so that the solution goes to infinity.
My issue is that I have been asked to plot the iterative solutions, but cant seem to understand how to handle the infinity, when storing the iterator and the respective solution values in a list inside the for loop, and when I call that list outside the loop. Below is the entire code and result pane with the error message.
import math
import numpy as np
import matplotlib.pyplot as plt
# Problem function #
def diffEq(x, y):
return (5*math.exp(y)) + (x)**2
# 4th order Runge - Kutta method #
def RK_fourth(x0,y0,x,h):
average_slope = 0
i = 0
x_axis = [] # For plotting #
y_axis = [] # For plotting #
print("Iterations of 4th Order Runge - Kutta Method")
for i in np.arange(x0,x,h, dtype=float):
k1 = diffEq(y0,i)
k2 = diffEq(y0 + (k1 * (h/2)), i + (h/2))
k3 = diffEq(y0 + (k2 * (h/2)), i + (h/2))
k4 = diffEq(y0 + (k3 * h), i + h)
average_slope = (1/6) * (k1 + (2* k2) + (2*k3) + k4)
print('average_slope:',average_slope)
y = y0 + (average_slope * h)
y0 = y
x_axis.append(i + h)
y_axis.append(y)
print("x:", i+h)
print("y:", y)
print(x_axis)
print(y_axis)
print(x_axis)
print(y_axis)
print("The solution using 4th order Runge - Kutta method is : ","%.4f"%y)
# main #
RK_fourth(0, 0 ,1, 0.1)
# Console #
Iterations of 4th Order Runge - Kutta Method
average_slope: 5.350250926709215
x: 0.1
y: 0.5350250926709216
[0.1]
[0.5350250926709216]
average_slope: 6.568756155777355
x: 0.2
y: 1.191900708248657
[0.1, 0.2]
[0.5350250926709216, 1.191900708248657]
average_slope: 9.107515760505036
x: 0.30000000000000004
y: 2.1026522842991606
[0.1, 0.2, 0.30000000000000004]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606]
average_slope: 14.993810753267782
x: 0.4
y: 3.602033359625939
[0.1, 0.2, 0.30000000000000004, 0.4]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939]
average_slope: 33.7279469045296
x: 0.5
y: 6.9748280500789
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789]
average_slope: 185.38472029098983
x: 0.6
y: 25.51330007917788
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788]
average_slope: 2568779.276837211
x: 0.7000000000000001
y: 256903.4409838003
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003]
average_slope: 1.4658111185680487e+68
x: 0.8
y: 1.4658111185680488e+67
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001, 0.8]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003, 1.4658111185680488e+67]
---------------------------------------------------------------------------
OverflowError Traceback (most recent call last)
/var/folders/bt/kqf88mw53h55m9mj35rkwt6h0000gn/T/ipykernel_18686/1130051705.py in <module>
6 #ImprovedEuler(0.2, 0.1 ,0.3, 0.1)
7 #print("\n")
----> 8 RK_fourth(0, 0 ,1, 0.1)
9 #print("\n")
10 #GaussJordan(12,3,-5,1,1,5,3,28,3,7,13,76) # Enter the coefficients #
/var/folders/bt/kqf88mw53h55m9mj35rkwt6h0000gn/T/ipykernel_18686/303333894.py in RK_fourth(x0, y0, x, h)
74 k1 = diffEq(y0,i)
75 k2 = diffEq(y0 + (k1 * (h/2)), i + (h/2))
---> 76 k3 = diffEq(y0 + (k2 * (h/2)), i + (h/2))
77 k4 = diffEq(y0 + (k3 * h), i + h)
78 average_slope = (1/6) * (k1 + (2* k2) + (2*k3) + k4)
/var/folders/bt/kqf88mw53h55m9mj35rkwt6h0000gn/T/ipykernel_18686/303333894.py in diffEq(x, y)
12 def diffEq(x, y):
13 #Function for GVF
---> 14 return (5*math.exp(y)) + (x)**2
15
16
OverflowError: (34, 'Result too large')
Note that I haven't called the plt.plot(x_axis, y_axis) yet. In order to verify that x_axis and y_axis is able to get all the values I called them once inside the loop and once outside the loop. The last print before the error message is from the print() statements inside the loop.
How to solve this issue, I am very much new to Python. I apologise if the issue is not clear to you and will be happy to clarify further.
Thank you in advance for helping out.
You can use numpy's np.exp instead of math.exp to handle infinity:
def diffEq(x, y):
return (5*np.exp(y)) + (x)**2
Result:
Iterations of 4th Order Runge - Kutta Method
average_slope: 5.350250926709215
x: 0.1
y: 0.5350250926709216
[0.1]
[0.5350250926709216]
average_slope: 6.568756155777355
x: 0.2
y: 1.191900708248657
[0.1, 0.2]
[0.5350250926709216, 1.191900708248657]
average_slope: 9.107515760505036
x: 0.30000000000000004
y: 2.1026522842991606
[0.1, 0.2, 0.30000000000000004]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606]
average_slope: 14.993810753267782
x: 0.4
y: 3.602033359625939
[0.1, 0.2, 0.30000000000000004, 0.4]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939]
average_slope: 33.7279469045296
x: 0.5
y: 6.9748280500789
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789]
average_slope: 185.38472029098983
x: 0.6
y: 25.51330007917788
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788]
average_slope: 2568779.276837211
x: 0.7000000000000001
y: 256903.4409838003
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003]
average_slope: 1.4658111185680487e+68
x: 0.8
y: 1.4658111185680488e+67
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001, 0.8]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003, 1.4658111185680488e+67]
average_slope: inf
x: 0.9
y: inf
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001, 0.8, 0.9]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003, 1.4658111185680488e+67, inf]
average_slope: inf
x: 1.0
y: inf
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001, 0.8, 0.9, 1.0]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003, 1.4658111185680488e+67, inf, inf]
[0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7000000000000001, 0.8, 0.9, 1.0]
[0.5350250926709216, 1.191900708248657, 2.1026522842991606, 3.602033359625939, 6.9748280500789, 25.51330007917788, 256903.4409838003, 1.4658111185680488e+67, inf, inf]
The solution using 4th order Runge - Kutta method is : inf
The program will output a runtime warning RuntimeWarning: overflow encountered in double_scalars.... You can suppress this warning by using a context manager:
import warnings
with warnings.catch_warnings():
warnings.filterwarnings("ignore", "overflow")
RK_fourth(0, 0 ,1, 0.1)