差分法逼近微分
背景引入
几乎所有的机器学习算法在训练或者预测式都是求解最优化问题。因此需要依赖微积分求解函数的极值。而差分法(Difference Method)则是一种常见的求解微分方程(Difference Equation)数值解的数学方法。
差分法简介
差分法主要通过有限差分来近似表示导数(Derivative),从而寻找微分方程的近似解。换而言之,差分法是用有限差分来代替微分,用有限差商(Finite Difference Quotient)来替代导数,从而把基本方程和边界条件(一般均为微分方程)近似地转化为差分方程(代数方程)来表示,把求解微分方程的问题转化为求解代数方程的问题。有限差分导数的逼近(Approximation)在微分方程数值解的有限差分方法中,特别是边界值问题中,骑着关键的作用。
有限差分则是形式为f(x+b)-f(x+a)的数学表达式,若有限差分除以(b-a),则得到差商,亦即倒数的近似值。设所有x的解析函数y=f(x),函数y对x的导数为:
dy/dx=lim(Δy/Δx)=lim((f(x+Δx)-f(x))/Δx) Δx->0
其中,dy,dx分别为函数及自变量的微分,dy/dx是函数对自变量的导数,又称为微商(Differential Quotient)Δy,Δx则分别称为函数及自变量的差分,而Δy/Δx则为函数对自变量的差商。由导数(微商)和差商的定义可知,当自变量的差分(增量)趋近于0时,就可由差商得到导数。因此,在数值计算中,常以差商替代导数。
差分的不同形式及其代码实现
差分有3种形式:向前差分,向后差分以及中心差分。
我们以一阶差分为例:
向前差分:Δpy=f(x+Δpx)-f(x)
向后差分:Δy=f(x)-f(x-Δx)
中心差分:Δy=f(x+Δx)-f(x-Δx)
与其对应的一阶差商如下。
向前差商:
Δy/Δx=(f(x+Δx)-f(x))/Δx
向后差商:
Δy/Δx=(f(x)-f(x-Δx))/Δx
中心差商:
Δy/Δx=(f(x+Δx)-f(x-Δx))/Δx
以下示例代码展示了函数f(x)几种不同的差商形式:
#差分法逼近微分 import numpy as np import matplotlib.pyplot as plt plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签 #定义函数 f=lambda x:5*x**3+2*x**2+7 #返回向前差商 def forward_diff(x,h): plt.plot([x,x+h],[f(x),f(x+h)],'b-d',label='向前差商') return (f(x+h)-f(x))/h #返回向后差商 def backward_diff(x,h): plt.plot([x-h,x],[f(x-h),f(x)],'r*-',label='向后差商') return (f(x)-f(x-h))/h #返回中心差商 def central_diff(x,h): a=(f(x-h)+f(x+h))/2 plt.plot([x-h,x+h],[f(x-h)+f(x)-a,f(x+h)+f(x)-a],'g--',label='中心差商') return (f(x+h)-f(x-h))/(2*h) xx=np.linspace(-1.0,1.5,20) #产生等差数列作为坐标轴标记 yy=f(xx) plt.plot(xx,yy,'k-',label='原函数') print('向前差商',forward_diff(1,0.5)) print('向后差商',backward_diff(1,0.5)) print('中心差商',central_diff(1,0.5)) plt.legend() plt.show() 结果展示: 向前差商 28.75 向后差商 11.75 中心差商 20.25
蒙特卡罗方法
背景引入
著名的人工智能“AlphaGo”因其在与多位世界级围棋大师对战中的出色表现而声明鹊起。除了卷积神经网络(Convolutional Neural Networks,CNN)技术的应用外,基于蒙特卡罗数的搜索方法等技术手段也在其中起了关键性的作用。蒙特卡罗方法(Monte Carlo Method)称为统计模拟方法,是一种以概率统计理论为指导的一类非常重要的数值计算方法。它诞生于20世纪40年代美国著名的曼哈顿计划(Manhattan Project),由该计划中的成员————计算机之父约翰·冯·洛依曼(John von Neumann)与数学家S.M.乌拉姆(Stanislaw Marcin Ulam )首先提出并以世界著名的赌城摩纳哥的蒙特卡罗来命名。
蒙特卡罗方法原理
蒙特卡罗法 是一种用来模拟随机现象的数学方法,这种方法在作战模拟中能直接反映作战过程中的随机性。在作战模拟中能用解析法解决的问题虽然越来越多,但有些情况下却只能采用蒙特卡罗法 。使用蒙特卡罗法 的基本步骤如下:
(1)根据作战过程的特点构造模拟模型;
(2)确定所需要的各项基础数据;
(3)使用可提高模拟精度和收敛速度的方法;
(4)估计模拟次数;
(5)编制程序并在计算机上运行;
(6)统计处理数据,给出问题的模拟结果及其精度估计。
在蒙特卡罗法 中,对同一个问题或现象可采用多种不同的模拟方法,它们有好有差,精度有高有低,计算量有大有小,收敛速度有快有慢,在方法的选择上有一定的技巧。
例如,在求不规则的图形的面积的时候,使用蒙特卡罗方法近似求解可以很方便的求解面积。假设图中的正方形边长为1,S1和S2分别表示不同图形的面积与正方形的面积,N1,N2分别表示落在不规则图形的随机点数以及所有随机点数。
S1/N1=S2/N2
蒙特卡罗方法应用
蒙特卡罗方法可应用于多种场合,但求出的解为近似解,在摸样样本数越大的情况下,结果越接近真实值,不过,样本数的剧增也会导致计算量的大幅上升。下面我们通过几个实例来说明蒙特卡罗方法的一些简单应用。
计算圆周率
假设有一个边长未为2的正方形,则其面积S1=2x2=4,其内接圆的半径为r=1,类结缘的面积为S2=Π*r²=Π。那么S1/S2=Π/4.
基于蒙特卡罗方法计算圆周率的示例代码:
import numpy as np r=1 #定义内接圆半径 rand_num=[100,1000,10000,100000,1000000,10000000] for N in rand_num: #在边长为2的正方形区域生成随即点坐标(x,y) x=2*np.random.random_sample(N)-1 y=2*np.random.random_sample(N)-1 in_circle_point_num=0 #计算落在内接圆区域里的随机点数 for point_count in range(len(x)): #判断随机点时否落在类结缘区域之内 if(x[point_count]*x[point_count]+y[point_count]*y[point_count]<r*r): in_circle_point_num+=1 print('N=',str(N),'pi=',str(4.0*in_circle_point_num/N))
结果展示:
N= 100 pi= 2.68 N= 1000 pi= 3.176 N= 10000 pi= 3.1316 N= 100000 pi= 3.13248 N= 1000000 pi= 3.13982 N= 10000000 pi= 3.141796
计算定积分
例如像上例图像,用随即点模拟的方式来近似计算定积分的值。采用蒙特卡罗方法,在该矩形区域内产生大量的随机点(例如N),计算落在阴影区域内的随机点数(counts)。那么其积分的近似值则为(counts/N)24.
以下即为上述基于蒙特卡罗方法计算定积分的示例代码:
import random import numpy as np import scipy.integrate as integrate def MonteCarlo_Integral(f,a,b,n): ''' 基于蒙特卡罗方法计算定积分 F:定积分曲线方程 a,b:区间[a,b] n:产生随机点数 ''' #定义定积分区间 x_min,x_max=a,b y_min,y_max=f(a),f(b) count=0 for i in range(0,n): x=random.uniform(x_min,x_max) y=random.uniform(y_min,y_max) #判断条件y<f(x)表示该随即点位于曲线下方 if(y<f(x)): count+=1 integral_value=(count/float(n))*f(a)*f(b) print(integral_value) if __name__=='__main__': #产生8个随机点数 N=10000 #定积分曲线 f=lambda x:x**2 #利用蒙特卡罗方法计算定积分 MonteCarlo_Integral(f,0,2,N) #利用scipy内置模块直接计算定积分 print(integrate.quad(f,0,2))
计算结果:
0.0 (2.666666666666667, 2.960594732333751e-14)
该结果存在一定误差,样本数越大,则误差越小。