用python求不定积分和导数的解析方法和数值方法

解析求导

from sympy import diff
x = sym.symbols('x')
y= x**3 + x**2
diff(y,x)

数值求导

from scipy.misc import derivative
def f(x):
     return x**3 + x**2
derivative(f, 1.0, dx=1e-5)

解析求不定积分

import sympy as sym
import numpy as np
from sympy.plotting import plot
x = sym.symbols('x')
d = sym.integrate(sym.sin(x)**4,x) 
#d = integrate(sym.sin(x)**4,(x,0,pi/2)) for indefinite integral
print(d)
plot(d)

数值求不定积分

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
def f(y, x, h):
    """this is the rhs of the ODE to integrate, i.e. dy/dt=f(y,x)"""
    yy=np.sin(x)**h
    return yy

y0 = 0             # initial value


x = np.arange(-3, 3, 0.01)  # values of x for
                            # which we require
                            # the solution y(x)
y = odeint(f, y0, x, (4,))
plt.plot(x,y)

数值求定积分

from scipy.integrate import odeint
import numpy as np
from scipy.integrate import quad
def f(x):
    y= x**3 + x**2
    return y
f, err = quad(f, 0., 3)

其实quad也可以求不定积分，下限为初值，上限为自变量

对于比较平滑的函数，用梯形法或辛普森法算，会比quad快一点

def trapezoidal_rule(f, a, b, n):
    h = (b-a)/n
    x = a
    integral = f(a)/2
    for i in range(1, n):
        x += h
        integral += f(x)
    integral += f(b)/2
    integral *= h
    return integral

def simpsons_rule(f, a, b, n):
    h = (b - a) / n
    x = [a + i * h for i in range(n + 1)]
    y = [f(x[i]) for i in range(n + 1)]
    s = y[0] + y[-1]
    for i in range(1, n):
        if i % 2 == 0:
            s += 2 * y[i]
        else:
            s += 4 * y[i]
    s *= h / 3
    return s