当前位置： 首页 > news >正文

Python数值计算（27）—— 数值微分

news 2024/10/22 5:20:13

1. 基本知识

微积分的核心是导数，在高数中，通常是利用极限的方式来定义导数的：

当\Delta x趋近于零时，就是该点处的导数：

\frac{\mathrm{d} y}{\mathrm{d} x}=\lim_{\Delta x \to 0}\frac{f(x_i+\Delta x)-f(x_i))}{\Delta x}

从几何角度看，曲线上点x_i处的导数，就是该点处切线的斜率，导数对于函数的变化率至关重要。

2. 具有O(h)误差的计算方法

根据前面的定义，虽然我们无法通过数值计算取得极限为零，但是当变化量足够小的时候，是可以通过差分的方法，估算处该点处的导数值的，通过对泰勒多项式的展开，以及差分选用的方式，分为了前向差分，中间差分和后向差分，以下是误差限为O(h)的1~4阶差分公式：

#differential method for functionimport numpy as npdef diff_h1(f,x,h=1e-2,order=1,method='central'):"""differentiation formula of order 1 to 4 using the forward, central, or backward method.error at O(h)Parameters----------f : functionThe function to be differentiated.x : floatThe point at which to differentiate the function.h : float, optionalThe step size for the finite difference method. The default is 1e-2.order : int, optionalThe order of the differentiation. The default is 1.method : str, optionalThe method to use for the finite difference method. The default is 'central'.Returns-------floatThe derivative of the function at the given point.Raises------ValueErrorIf the order is greater than 4or if parameter h is less than 1e-3(due to floating point error).or if the method is not 'forward', 'central', or 'backward'.""" if order==0:return f(x)if h<1e-3:raise ValueError("h must be greater than 1e-3,due to floating point error")if order>=5:raise ValueError("Order must be less than 5")if method not in ['forward','central','backward']:raise ValueError("Method must be 'forward', 'central', or 'backward'")if method == 'forward':if order==1:return (f(x+h)-f(x))/helif order==2:return (f(x+2*h)-2*f(x+h)+f(x))/(h**2)elif order==3:return (f(x+3*h)-3*f(x+2*h)+3*f(x+h)-f(x))/(h**3)elif order==4:return (f(x+4*h)-4*f(x+3*h)+6*f(x+2*h)-4*f(x+h)+f(x))/(h**4)if method == 'central':if order==1:return (f(x+h)-f(x-h))/(2*h)elif order==2:return (f(x+h)-2*f(x)+f(x-h))/(h**2)elif order==3:return (f(x+2*h)-2*f(x+h)+2*f(x-h)-f(x-2*h))/(2*h**3)elif order==4:return (f(x+2*h)-4*f(x+h)+6*f(x)-4*f(x-h)+f(x-2*h))/(h**4)if method == 'backward':if order==1:return (f(x)-f(x-h))/helif order==2:return (f(x)-2*f(x-h)+f(x-2*h))/(h**2)elif order==3:return (f(x)-3*f(x-h)+3*f(x-2*h)-f(x-3*h))/(h**3)elif order==4:return (f(x)-4*f(x-h)+6*f(x-2*h)-4*f(x-3*h)+f(x-4*h))/h/h/h/hraise ValueError("Unknown method or order")

注意：上述公式中，计算差分的公式中，h的值不能小于0.0001，因为很有可能此时在高阶（通常是4阶）差分的时候，由于浮点数的精度带来的影响，计算结果会出现较大偏差，这并不是方法本身的问题，而是由于浮点数计算精度导致的。

3. 具有O(h^2)误差的计算方法

通过上述各公式，可知在需要计算1阶差分的时候，需要计算两个函数值，计算4阶差分时，就需要用到五个点处的函数值，但是如果通过引入更多的点，实际上是可以构造出误差精度为O(h^2)的1~4阶差分公式，如下所示：

def diff_h2(f,x,h=1e-2,order=1,method='central'):"""differentiation formula of order 1 to 4 using the forward, central, or backward method.error at O(h^2)Parameters----------f : functionThe function to be differentiated.x : floatThe point at which to differentiate the function.h : float, optionalThe step size for the finite difference method. The default is 1e-2.order : int, optionalThe order of the differentiation. The default is 1.method : str, optionalThe method to use for the finite difference method. The default is 'central'.Returns-------floatThe derivative of the function at the given point.Raises------ValueErrorIf the order is greater than 4 or if parameter h is less than 1e-3(due to floating point error).or if the method is not 'forward', 'central', or 'backward'.""" if order==0:return f(x)if h<1e-3:raise ValueError("h must be greater than 1e-3,due to floating point error")if order>=5:raise ValueError("Order must be less than 5")if method not in ['forward','central','backward']:raise ValueError("Method must be 'forward', 'central', or 'backward'")if method == 'forward':if order==1:return (-f(x+2*h)+4*f(x+h)-3*f(x))/(2*h)elif order==2:return -(f(x+3*h)+4*f(x+2*h)-5*f(x+h)+2*f(x))/(h**2)elif order==3:return (-3*f(x+4*h)+14*f(x+3*h)-24*f(x+2*h)+18*f(x+h)-5*f(x))/(2*h**3)elif order==4:return (-2*f(x+5*h)+11*f(x+4*h)-24*f(x+3*h)+26*f(x+2*h)-14*f(x+h)+3*f(x))/(h**4)if method == 'central':if order==1:return (-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h)elif order==2:return (-f(x+2*h)+16*f(x+h)-30*f(x)+16*f(x-h)-f(x-2*h))/(12*h**2)elif order==3:return (-f(x+3*h)+8*f(x+2*h)-13*f(x+h)+13*f(x-h)-8*f(x-2*h)+f(x-3*h))/(8*h**3)elif order==4:return (-f(x+3*h)+12*f(x+2*h)-39*f(x+h)+56*f(x)-39*f(x-h)+12*f(x-2*h)-f(x-3*h))/(6*h**4)if method == 'backward':if order==1:return (3*f(x)-4*f(x-h)+f(x-2*h))/(2*h)elif order==2:return (2*f(x)-5*f(x-h)+4*f(x-2*h)-f(x-3*h))/(h**2)elif order==3:return (5*f(x)-18*f(x-h)+24*f(x-2*h)-14*f(x-3*h)+3*f(x-4*h))/(2*h**3)elif order==4:return (3*f(x)-14*f(x-h)+26*f(x-2*h)-24*f(x-3*h)+11*f(x-4*h)-2*f(x-5*h))/h/h/h/hraise ValueError("Unknown method or order")

4. 测试

测试如下：

假设有原函数：

f(x)=e^{-x}sinx

其4阶导函数为：

f(x)=-4e^{-x}sinx

在此处需要提出注意，就是在很多时候我们并不一定能够获取到实际的解析函数关系式，但是又需要计算采样点附近的导数值，以上给出一个具体的具有解析形式的函数，主要是为了演示差分逼近导数的算法。

以下测试代码输出真实导数值，以及通过差分方式给出的近似结果如下：

x=0.5
print(df4(x))
h=1e-3
print(diff_h2(f,x,h,4,'central'))
# output:
#-1.1631451528507675
#-1.1632546777680846

可以看到两者已经非常接近，其差为0.0001095249173171，似乎并没有达到号称的O(h^2)误差范围，这个问题还是由于分母太小的缘故。

如果选择h=0.1：

x=0.5
print(df4(x))
h=1e-1
print(diff_h2(f,x,h,4,'central'))
# output:
#-1.1631451528507675
#-1.1631586495957265

这次两者的差为0.000013496744959，可以看到实际上已经达到了0.00001的误差水平，不仅比前面的误差还小，甚至还比O(h^2)误差范围还小很多。

所以，在数值计算的时候，虽然公式可以保证，较小的步长可以产生更加精确的结果，但是由于在实际计算时，由于浮点数计算的限制，以及分母太小，容易造成分子的误差被放大，最终结果是更小的步长并不一定会产生更小的误差。


http://www.mrgr.cn/news/55733.html

相关文章：

  • 【Eclipse系列】The word is not correctly spelled问题解决
  • 解决Git合并冲突：掌握版本控制的精髓
  • C++面向对象编程学习
  • Qt编写的modbus模拟器/支持网络和串口以及websocket/支持网络rtu
  • 使用finalshell远程ssh连接ubuntu
  • Qt中使用线程之QThread
  • 基于Springboot在线视频网站的设计与实现
  • 心觉：突破自己
  • 51单片机快速入门之 IIC I2C通信
  • UML之用例图详解
  • 【ShuQiHere】深入了解逻辑门与晶体管数量：CMOS技术详解
  • 毕业设计选题：基于Hadoop的热点新闻分析系统的设计与实现
  • js构造函数和原型对象，ES6中的class，四种继承方式
  • Python Flask 数据库开发
  • 提示词高级阶段学习day3.1
  • 目前最新 Reflector V11.1.0.2067版本 .NET 反编译软件
  • 【C++】拆分详解 - stack和queue
  • 03_深入理解Linux：系统组成、内核版本及文件系统详解
  • 【MySQL】索引和事务
  • JAVA继承
  • 时间数据可视化基础实验——Python实现
  • 【付费】Ambari集成Dolphin实战-002-bigtop下编译dolphin——下
  • 简述 C# 二维数据集合 List 的创建、遍历、修改、输出
  • 3. IoC 与DI
  • 数据流风格
  • 改变函数调用上下文：apply与call方法详解及实例