课程概要：
　　1、非线性方程组
　　2、数值积分
　　3、常微分方程组

---

1、非线性方程组

'''
    5 * x1 - 25 = 0
    5 * x0 * x0 - x1 * x2 = 0
    x2 * x0 -27 = 0

scipy.optimize fsolve(func, x)          
#   所使用的scipy中的库optimize以及方法fsolve
#   func 是自己构造的函数
'''

from scipy.optimize import fsolve

def func(x):
    x0, x1, x2 = x.tolist()
    return [
            5 * x1 - 25,
            5 * x0 * x0 - x1 * x2,
            x2 * x0 -27
           ]

r = fsolve(func,[1,1,1])
print r

>>> 
[ 3.  5.  9.]

'''
    5 * x1  + 3 = 0
    5 * x0 * x0 - 2sin(x1 * x2) = 0
    x1 * x2 -1.5 = 0
'''

from scipy.optimize import fsolve
from math import sin

def func(x):
    x0, x1, x2 = x.tolist()
    return [
                5 * x1  + 3,
                5 * x0 * x0 - 2 * sin(x1 * x2) ,
                x1 * x2 -1.5
                ]

r = fsolve(func,[1,1,1])
print r

>>> 
[-0.63166288  -0.6   -2.5 ]

'''
    5 * x1  + 3 = 0
    5 * x0 * x0 - 2sin(x1 * x2) = 0
    x1 * x2 -1.5 = 0
'''

from scipy.optimize import fsolve
from math import sin, cos

def func(x):
    x0, x1, x2 = x.tolist()
    return [
                5 * x1  + 3,
                5 * x0 * x0 - 2 * sin(x1 * x2) ,
                x1 * x2 -1.5
                ]

def j(x):
    x0, x1, x2 = x.tolist()
    return [
          [0, 5, 0],        #   分别对x0, x1, x2 求偏导
          [10* x0, -2 * x2 * cos(x1 * x2), -2 * x1 * cos(x1 * x2)],
          [0, x2, x1]
           ]

r = fsolve(func, [1,1,1], fprime=j)     #   使用求导矩阵传入时可以提高4倍效率
print r

>>>
 [-0.63166288   -0.6    -2.5]

---

2、数值积分（定积分）

#   求半圆面积（定义求）
'''
    pi * 1 *1 / 2

    x * x + y * y = 1
    y = (1 - x * x) ** 0.5
'''

import numpy as np

def func(x):
    return (1 - x * x) ** 0.5

N = 10000
x = np.linspace(-1, 1, N)
dx = 2.0 / N
y = func(x)

m = dx * np.sum(y[:-1] + y[1:])
print m / 2                     #   m是整圆的面积

>>> 
1.570637584

#   使用scipy的库进行计算
import numpy as np
from scipy import integrate

def func(x):
return (1 - x * x) ** 0.5           #   x^2 + y^2 = 1
    
    p, err = integrate.quad(func, -1, 1)        #   x 给定的区间
print p

>>> 
1.57079632679

#   使用scipy的库进行计算
import numpy as np
from scipy import integrate


def func(x):
return (4 - x * x) ** 0.5       #   x^2 + y^2 = 4
    
    p, err = integrate.quad(func, -2, 2)
print p

>>> 
6.28318530718

#   dblquad：双重积分        trlquad：三重积分

'''
    x * x+ y * y + z * z =1
    z = (1 - x * x - y * y)**0.5
'''

from scipy import integrate

def func(x, y):
    return (1 - x * x - y * y)**0.5

def func2(x):
    return (1 - x * x)**0.5

m = integrate.dblquad(func, -1, 1,          #   x 的区间
                       lambda x: - func2(x),
                       lambda x: func2(x))      #   y 的积分区间

print m

>>> 
(2.0943951023931984, 2.3252456653390912e-14)

---

3、常微分方程组

dx/dt=σ(y-x)
dy/dt=x(ρ-z)-y
dz/dt=xy-βz
#   σ = p , ρ = r , β= b

from scipy.integrate import odeint
import numpy as np

def lorenz(w, t, p, r, b):
    x, y, z = w
    return np.array([p * (y-x), x * (r-z) - y, x * y - b * z])

t = np.arange(0, 30, 0.01)
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))        #   初始值，即w
##  track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))    #   args 为(p, r, b)

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:, 0], track1[:, 1], track1[:, 2])
plt.show()

#   改变轨迹值
ax.plot(track2[:, 0], track2[:, 1], track2[:, 2]) 
#   只有细微的差别，当改为1.10时差距就会很明显，但第二个图与视频中绘出的图有很大区别