Lab 9 – Taylor method
from numpy import array,zeros,exp
def taylor ( deriv ,x ,y , xStop , h ) :
X = [ ]
Y = [ ]
X . append ( x )
Y . append ( y )
while x < xStop :
D = deriv (x , y )
H = 1.0
for j in range ( 3 ) :
H = H * h / ( j + 1 )
y = y + D [ j ] * H # H = h ^ j / j !
x = x + h
X . append ( x ) # Append results to
Y . append ( y ) # lists X and Y
return array ( X ) , array ( Y ) # Convert lists into arrays
def deriv (x , y ) :
D = zeros (( 4 , 1 ) )
D [ 0 ] = [ 2 * y [ 0 ] + 3 * exp ( x ) ]
D [ 1 ] = [ 4 * y [ 0 ] + 9 * exp ( x ) ]
D [ 2 ] = [ 8 * y [ 0 ] + 21 * exp ( x ) ]
D [ 3 ] = [ 16 * y [ 0 ] + 45 * exp ( x ) ]
return D
x = 0.0 # Initial value of x
xStop = 0.3 # last value
y = array ( [0.0 ] ) # Initial values of y
from numpy import array, zeros, exp
def taylor ( deriv ,x ,y , xStop , h ) :
X = [ ]
Y = [ ]
X . append ( x )
Y . append ( y )
while x < xStop : # Loop over integration steps
D = deriv (x , y ) # Derivatives of y
H = 1.0
for j in range ( 3 ) : # Build Taylor series
H = H * h / ( j + 1 )
y = y + D [ j ] * H # H = h ^ j / j !
x = x + h
X . append ( x ) # Append results to
Y . append ( y ) # lists X and Y
return array ( X ) , array ( Y ) # Convert lists into arrays
# deriv = user - supplied function that returns derivatives in the 4 x
n
# array
def deriv (x , y ) :
D = zeros (( 4 , 1 ) )
D [ 0 ] = [ 2 * y [ 0 ] + 3 * exp ( x ) ]
D [ 1 ] = [ 4 * y [ 0 ] + 9 * exp ( x ) ]
D [ 2 ] = [ 8 * y [ 0 ] + 21 * exp ( x ) ]
D [ 3 ] = [ 16 * y [ 0 ] + 45 * exp ( x ) ]
return D
x = 0.0 # Initial value of x
xStop = 0.3 # last value
y = array ([ 0.0 ] ) # Initial values of y
h = 0.1 # Step size
X , Y = taylor ( deriv ,x ,y , xStop , h )
print ( " The required values are : at x = %.2f , y = %.5f , x = %.2f ,
y = %.5f , x = %.2f , y = %.5f , x = %.2f , y = %.5f " % ( X [ 0 ] ,Y [
0 ] ,X [ 1 ] ,Y [ 1 ] ,X [ 2 ] ,Y [ 2 ] ,X [ 3 ] ,Y [ 3 ] ) )
1. The ODE:
Although not explicitly stated in the code, the ODE being solved can be inferred from
the deriv function:
library_add
content_copy
def deriv(x, y):
D = zeros((4, 1))
D[0] = [2 * y[0] + 3 * exp(x)] # First derivative
D[1] = [4 * y[0] + 9 * exp(x)] # Second derivative
D[2] = [8 * y[0] + 21 * exp(x)] # Third derivative
D[3] = [16 * y[0] + 45 * exp(x)] # Fourth derivative
return D
Use code with caution
This function calculates the first four derivatives of the unknown function y(x). From
the first derivative, we can see that the ODE is likely of the form:
library_add
content_copy
dy/dx = 2*y + 3*exp(x)
Use code with caution
2. The Taylor Method:
The Taylor method uses the Taylor series expansion of the solution function y(x) to
approximate its value at a given point. The Taylor series expansion around a
point x is:
library_add
content_copy
y(x + h) = y(x) + h*y'(x) + (h^2/2!)*y''(x) + (h^3/3!)*y'''(x) + ...
Use code with caution
where:
h is the step size
y'(x), y''(x), y'''(x), ... are the derivatives of y(x) evaluated at x.
3. Implementation in the Code:
The taylor function implements the Taylor method.
It initializes two lists, X and Y, to store the values of x and y at each step.
It iterates until x reaches xStop, performing the following steps in each
iteration:
o Calculates the derivatives of y using the deriv function.
o Builds the Taylor series expansion up to the third-order term,
accumulating the terms in y.
o Updates x by adding the step size h.
o Appends the new values of x and y to the lists X and Y.
Finally, it converts the lists X and Y into NumPy arrays and returns them.
4. The Specific Problem:
The code aims to solve the ODE (likely dy/dx = 2*y + 3*exp(x)) with the initial
condition y(0) = 0 over the interval x = 0 to x = 0.3 using a step size of h = 0.1.
The output displays the values of y at x = 0.0, x = 0.1, x = 0.2, and x = 0.3.
In essence, the code is using a third-order Taylor series approximation to numerically
solve the given ODE and find the values of the solution at specific points.
Lab 9 – Eulers method
import numpy as np
import matplotlib.pyplot as plt
# Define parameters
f = lambda x, y:np.exp(-x) # ODE
h = 0.2 # Step size
y0 = -1 # Initial Condition
n = 3
# Explicit Euler Method
x = np.zeros(n + 1) # Initialize x as an array
y = np.zeros(n + 1) # Initialize y as an array
y[0] = y0
x[0] = 0
for i in range(0, n):
x[i + 1] = x[i] + h
y[i + 1] = y[i] + h * f(x[i], y[i])
print("The required values are at x = % 0.2f , y = % 0.5f , x = %
0.2f , y = % 0.5f,x = % 0.2f , y = % 0.5f , x = % 0.2f , y = % 0.5f " %
(x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]))
#print("The required values are at x = % 0.2f , y = % 0.5f , x = %
0.2f ,
#y = % 0.5f,x = % 0.2f , y = % 0.5f , x = % 0.2f ,
#y = % 0.5f " % (x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]))
print(" \ n \ n ")
plt.plot(x, y, 'bo--', label=' Approximate ')
plt.plot(x -np.exp(-x), 'g*--', label='Exact')
plt.title(" Approximate and Exact Solution ")
plt.xlabel('x')
plt.ylabel('f( x )')
plt.grid()
plt.legend(loc='best')
plt.show()
1. The Ordinary Differential Equation (ODE):
The code aims to solve an ODE represented by the function:
library_add
content_copy
f = lambda x, y: np.exp(-x)
Use code with caution
This defines a first-order ODE where the derivative of y with respect to x (dy/dx) is
equal to exp(-x). In mathematical notation:
library_add
content_copy
dy/dx = exp(-x)
Use code with caution
2. Explicit Euler Method:
The code uses the Explicit Euler method to approximate the solution to this ODE.
This method is a numerical technique for solving differential equations. Here's how it
works:
Discretization: The method divides the x-axis into discrete steps of size h.
Approximation: It uses the following formula to approximate the value
of y at the next step (x + h):
library_add
content_copy
y(x + h) ≈ y(x) + h * f(x, y(x))
Use code with caution
where f(x, y(x)) is the value of the derivative at the current step.
Iteration: This process is repeated for multiple steps to obtain an
approximate solution curve.
3. The Exact Solution:
The exact solution to the ODE dy/dx = exp(-x) can be found by integrating both
sides:
library_add
content_copy
∫ dy = ∫ exp(-x) dx
Use code with caution
This gives us:
library_add
content_copy
y = -exp(-x) + C
Use code with caution
where C is the constant of integration. Using the initial condition y(0) = -1, we can
solve for C:
library_add
content_copy
-1 = -exp(0) + C
-1 = -1 + C
C = 0
Use code with caution
Therefore, the exact solution is:
library_add
content_copy
y = -exp(-x)
Use code with caution
4. Exponential Decay:
The function -exp(-x) is an example of exponential decay. Here's why:
Base of the exponent: The base of the exponent is e (Euler's number),
which is approximately 2.718.
Negative exponent: The exponent is -x. As x increases, the value of -
x becomes more negative.
Decaying behavior: Since e is raised to a negative and increasingly
negative power, the overall value of the function decreases as x increases.
Asymptote: The function approaches 0 as x approaches infinity but never
actually reaches 0. This is why the graph has a horizontal asymptote at y = 0.
