import math import time # import time to use for performance analysis import numpy as np # import numpy for array space import matplotlib.pyplot as plt # import matplotlib for graphing functions from scipy.integrate import odeint # import scipy to use the ordinary differencial equation integral function # Ordinary Differential Equation def dy_dx(y, x): # takes inputs y and x return 3*x + 2*y # returns the function y/(e^x - 1) # Runge-Kutta Forumla based off the formula used in class def rungeKutta(y, x, h): # takes y, x, and h inputs k1 = dy_dx(y, x) # solves k1 using the differential equation function k2 = dy_dx(y + ((0.5 * h) * k1), x + 0.5 * h) # solves k2 based on the answer from k1 using the differential equation function k3 = dy_dx(y + ((0.5 * h) * k2), x + 0.5 * h) # solves k3 based on the answer from k2 using the differential equation function k4 = dy_dx(y + (h * k3), x + h) # solves k4 based on the answer from k3 using the differential equation function t4 = (1.0 / 6.0) * (k1 + (2 * k2) + (2 * k3) + k4) # solves for t4 by taking a 6th of k1 + 2*k2 + 2*k3 + k4 y = y + (t4 * h) # solves for y by taking the initial y value and adding it to t4*h return y # returns y # initial variable values y = 4 # initial y value of 5 as outlined in the assignment x = 0 # initial x value of 1 as outlined in the assignment h = 0.02 # initial h value of 0.02 as outlined in the assignment n = 2000 # initial n value of 2000 chosen between 1000 or 2000 as outlined in the assignment # Runge-Kutta x and y arrays xsr = [] # x-space runge-kutta array to store the values of x for the runge-kutta function ysr = [] # y-space runge-kutta array to store the values of y for the runge-kutta function # ODEint x and y arrays, solution, and time analysis tso = time.time() # time start for ODEint function solution xso = np.linspace(0.0, 0.2 , 30) # x-space ODEint from 1 to n*h (40) plus 1 with a step size of n yso = odeint(dy_dx, y, xso) # y-space ODEint useing the odeint function from scicpy to find the y-space teo = time.time() # time end fore ODEint function solution tto = teo - tso # total time the ODEint function to solve # graphing ODEint plt.title("ODE Function Analysis") # set the title of the graph plt.xlabel("x") # set the x label on the graph plt.ylabel("y") # set the y label on the graph plt.plot(xso, yso, 'g-', label = "ODEint", linewidth = 2) # set the ODE line to be red and label it plt.legend() # shows the legend on the graph plt.savefig('Exact Solution1.png') # displays the graph # Runge-Kutta solution and time analysis tsr = time.time() # time start for runge-kutta function solution while (x <= 0.2): # for loop to run the runge-kutta function n number of times xsr.append(x) # append the x value to the x-space runge-kutta array ysr.append(y) # append the y value to the y-space runge-kutta array y = rungeKutta(y, x, h) # update the y value using the rungeKutta function x += h # update the x value by moving one step forward (0.02) ter = time.time() # time end for runge-kutta function solution ttr = ter - tsr # total time the runge-kutta function to solve td = ttr - tto # time difference between ODEint and runge-kutta function # graphing runge-kutta plt.title("Runge-Kutta Function Analysis") # set the title of the graph plt.xlabel("x") # set the x label on the graph plt.ylabel("y") # set the y label on the graph plt.plot(xsr, ysr, 'b--', label = "Runge Kutta") # set the runge-kutta to be blue and label it plt.legend() # shows the legend on the graph plt.savefig('Runge_Kutta1.png') # displays the graph # solutions print("\nODEint Solution: ", yso[-1]) # ODEint function solution print("Runge-Kutta Solution: ", ysr[-1]) # Runge-Kutta function solution # Print statement for time difference print("\nODEint Time: ", tto) # print the ODEint time print("Runge Kutta Time: ", ttr) # print the runge-kutta time print("ODEint is ", td, " seconds faster\n\n") # print the difference between ODEint and runge-kutta # error calculation error = 0 # initial error value of 0 errorRange = [] # array to store error over xn errorSpace = np.linspace(0.0,0.2, 30) # error space for error analysis for i in range(len(ysr)): # for loop to run through every x values error += (np.abs(ysr[i] - yso[i])/yso[i]) * 100 # sum all the error values using the percentage error formula errorRange.append((np.abs(ysr[i] - yso[i])/yso[i]) * 100) print("Percent Error at x =", i, ":", (np.abs(ysr[i] - yso[i])/yso[i]) * 100) # print error at each x value print("\nAverage Error Percent:", error/((int)(n * h) + 1), "\n") # print the total error divided by the total number of x values # graphing error #plt.title("Error Analysis") # set the title of the graph #plt.xlabel("xn") # set the x label on the graph #plt.ylabel("error") # set the y label on the graph #plt.plot(errorSpace, errorRange, label = "Error over Xn") # create the line and label it #plt.legend() # shows the legend on the graph #plt.show() # displays the graph # # graphing both functions plt.title("Runge-Kutta and ODE Function Analysis") # set the title of the graph plt.xlabel("x") # set the x label on the graph plt.ylabel("y") # set the y label on the graph plt.plot(xso, yso, 'r-', label = "ODEint", linewidth = 2) # set the ODE line to be red and label it plt.plot(xsr, ysr, 'bo', label = "Runge Kutta") # set the runge-kutta to be blue and label it plt.legend() # shows the legend on the graph plt.savefig('Comparison1.png')