Runge Kutta program outputting incorrect graphs

Question

I'm creating a program to compare two graphs, using the Runge Kutta method to solve an ODE. Ideally the resulting plots would be directly on top of each other, but the Runge Kutta function is outputting an incorrect plot. Am I incorrectly populating the array, or is it a problem with calling the new values?

Here is the program I have:

#Import correct libraries and extensions
import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
import warnings
def fxn():
    warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()

#Define conditions and store values
h=0.02 #Given values for ODE
x0=1
y0=5
xpoints=[x0] #array for storing X values
ypoints=[y0] #Array for storing Y values

#State equations and define functions
def dy_dx(y,x):
    return x / (np.exp(x) - 1) #Provided Equation
def RungeKuttaFehlberg(x,y):
    return x / (np.exp(x) - 1)
#Calculates k1-k4, x and y solutions
def RKFAlg(x0,y0,h):
    k1 = RungeKuttaFehlberg(x0,y0)
    k2 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k1))
    k3 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k2))
    k4 = RungeKuttaFehlberg(x0+h,y0+(h*k3))
    y1 = y0+(h/6)*(k1+(2*k2)+(2*k3)+k4)
    x1 = x0+h
    x1 = round(x1,2)
    print("Point (x(n+1),y(n+1)) =",(x1,y1))
    return((x1,y1)) #Returns as ordered pair

#Define range for number of calculations
for i in range(2000):
    print(f"Y{i+1}".format(i)) #Solution value format
    x0,y0 = RKFAlg(x0,y0,h) #Calls RKF Function
    xpoints.append(x0) #Saves values into array
    ypoints.append(y0)
    y0 = 1
    ODEy1 = odeint(dy_dx,y0,xpoints)

#Runge-Kutta Graph
plt.plot(xpoints,ypoints,'b:',linewidth = 1) #command to plot lines using various colors and widths
plt.suptitle("RKF Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()

#ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=1)
plt.suptitle("ODE Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()

#Function for plotting RKF and ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=2,label="ODE")
plt.plot(xpoints,ypoints,'b:',linewidth=3,label="Runge-Kutta")
plt.suptitle("ODE and RKF Comparison")
plt.legend(bbox_to_anchor=(.8,1),loc=0,borderaxespad=0)
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

#Function for plotting the difference graph
diff = [] #array to store difference
for i in range(len(xpoints)):
    diff.append(ypoints[i]-ODEy1[i])
plt.plot(xpoints,diff)
plt.suptitle("Difference")
plt.xlabel("x Points")
plt.ylabel("RKF and ODE diff.")
plt.show()

You have some strange names, there is no Fehlberg in sight, but that is not an error. In the main loop you have an indentation error, resetting `y0` should not be in the loop. — Lutz Lehmann, Oct 01 '22 at 06:19

score 0 · Answer 1 · answered Dec 15 '22 at 15:42

In python indentation defines the structure of a program, what the command block is inside a loop and what comes after the loop.

Your aim is to iterate through the time interval with steps h building the arrays for x and y. After this loop you want to call odeint once to get reference values over the x array.

Thus change this part to

#Define range for number of calculations
x,y = x0,y0  # for systems use y0.copy()
for i in range(2000):
    print(f"Y{i+1}".format(i)) #Solution value format
    x,y = RKFAlg(x,y,h) #Calls RKF Function
    xpoints.append(x) #Saves values into array
    ypoints.append(y)

ODEy1 = odeint(dy_dx,ypoints[0],xpoints)

For consistency it would be better to keep the interfaces of the integrators close,

def RK4Alg(func,x0,y0,h):
    k1 = func(x0,y0)
   ...

Then you can pass the same rhs function dy_dx to both methods, and can remove the mis-named RungeKuttaFehlberg.

For a sensible reference solution one should set the error tolerances to values that are decisively lower than the expected errors of the method to test.

opts = {"atol":1e-8, "rtol":1e-11}
ODEy1 = odeint(dy_dx,y0,xpoints,**opts)

Runge Kutta program outputting incorrect graphs

1 Answers1