How to resolve integration function not integrating correctly?

Question

I am trying to build a few simple operations such as a derivative and integral function to operate on lambda functions because sympy and scipy were struggling to integrate some things that I was passing to them.

The derivative function does not give me any issues and looks to return the derivative of the input function when plotted, but the integral function does not return the same, and does not plot the correct integral of the input.

import matplotlib.pyplot as plt
import numpy as np
from phys_func import func

sr = [-10,10]
x = np.linspace(sr[0],sr[1], 100)

F = lambda x: x**2
f = func.I(F,x)

plt.plot(x,F(x), label = 'F(x)')
plt.plot(x,f(x), label = 'f(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

The integration function that does not work:

def I(F,x):
        dx = (x[len(x)-1] - x[0])/len(x)
        return lambda x : 0.5*( F(x+dx) + F(x) )*dx

The derivative function that works:

def d(f,x):
        dx = (x[len(x)-1]-x[0])/len(x)
        return lambda x: (f(x+dx)-f(x))/dx

Can anyone lend me a hand please?

Answer 1

You cannot find the anti derivative of a function numerically with out knowing the value of the anti derivative at a single point. Suppose if you fix the value of the antiderivative function value at x =a to be 0 (and given function is continuous from [a,x]) , then we can use definite integrals. For this function, let us take a=0 (i.e 0 is a root of anti derivative function), so you can do a definite integral from [0,x]. Also, your integration function is wrong. You need to sum all the 0.5*( F(x+dx) + F(x) )*dx elements from 0 to x to get the definite integral. You can modify I(f,x) as follows

def I(F1): # N is number of intervals
    return lambda x, N: np.sum( 0.5*( F1(np.linspace(0,x,num=N)+ (x)/N ) + F1(np.linspace(0,x,num=N)))*(x)/N)

In [1]: import numpy as np

In [2]: def I(F1): # N is number of intervals
    ...:     return lambda x, N: np.sum( 0.5*( F1(np.linspace(0,x,num=N)+ (x)/N
    ...: ) + F1(np.linspace(0,x,num=N)))*(x)/N)
    ...:

In [3]: F = lambda x: x**2

In [4]: x_ran = np.linspace(-10,10, 100)

In [5]: y = I(F)

In [6]: y_ran = []

In [7]: for i in x_ran:
    ...:     y_ran.append(y(i,100))

In [8]: plt.plot(x_ran,y_ran)

In [9]: plt.show()