'''
Use the Taylor Series method of order 4 to compute $\int_0^2 e^-s^2 ds
by solving the Initial Value problem on the interval t \in [0,2]
Given dx/dt = e^-t^2
and   x(0) = 0
Check from the table of error function x(2) = (pi^(1/2)/2). erf(2) = 0.8820813907
We need to find the derivative till the fourth order of the Given function, and add them to a Taylor
Series to form the solution
'''

# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them



# we will begin by defining the function 
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return math.cos(t) - math.sin(x) + t**2
 

# now we differentiate the same function and find the 
# next orders of derivatives
def f2(t,x):
    return -math.sin(t) - math.cos(x)*f1(t,x) + 2*t

def f3(t,x):
    return -math.cos(t) + 2 - f2(t,x) * math.cos(x) + math.sin(x)*f1(t,x)**2

def f4(t,x):
    return math.sin(t) - (f2(t,x)*f1(t,x)*(-math.sin(x)) + math.cos(x)* f3(t,x)) + (math.cos(x)* f1(t,x)**3 + math.sin(x) * 2* f2(t,x) )


# also the function that finds the Error

def error(dt,m):
    return ((dt**m)/factorial(m+1)) * (xx[len(xx)-1] - xx[len(xx)-2])



# and the exact solution to the problem(for reference)
#def exact(t):
#    return math.sqrt(np.pi)*math.erf(t)*0.5



# Now we write down the values known to us
m = 4                       # order parameter
dt = 0.1                    # step Size
t,x = -1 ,3                    # initial values of both
tf = 2                      # final point to reach
xe = 0
# initialised arrays for graph plotting
tt = []
xx = []
xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x) + np.power(dt,4)/24*f4(t,x)       # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1





# Plotting and Visualization

#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'o',linewidth=2,color='red', label ='Taylor Method Order = 4')
#plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method_2.png')



# checking for the value at t = 2
# writing all the values 

#print("The Value that x(t) takes for different t:\n", xx)

# printing the error

#print("The Error Value is:\n ",error(dt,m))

#  difference between exact value and Taylor approximated

#print("The Difference between the Exact solution at t = 2 and Taylor approximated is:\n", math.erf(2)-xx[len(xx)-1])