【数值计算方法】数值积分微分-python实现-p3

原文链接：https://www.cnblogs.com/aksoam/p/18332123
更多精彩，关注博客园主页，不断学习！不断进步！

我的主页

csdn很少看私信，有事请b站私信

博客园主页-发文字笔记-常用

有限元鹰的主页 内容：

1. ABAQUS数值模拟相关
2. Python科学计算
3. 开源框架，编程学习笔记

哔哩哔哩主页-发视频-常用

FE-有限元鹰的个人空间 内容：

1. 模拟案例
2. 网格划分
3. 游戏视频，及其他搬运视频

6. 多重积分

本节考虑多重积分的两种计算方法:

• 第一种方法可看成前面几节关于单重积分的推广 ,计算重数较小的积分时有效 .
• 第二种方法称为蒙特卡罗 (Monte Carlo) 方法 , 是一种随机算法 , 处理高维积分时特别有效

6.1 从一重积分推广到二重积分

以二重积分为例,考虑:

$$I={\iint }_{\Omega }f(x,y)\mathrm{d}x\mathrm{d}y,$$

其中,$\Omega =a\le x\le b,{\phi }_1(x)\le y\le {\phi }_2(x)$

根据二重积分性质:

$$I={\iint }_{\Omega }f(x,y)\mathrm{d}x\mathrm{d}y={\int }_a^b\mathrm{d}x{\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)\mathrm{d}y.$$

令:

$$F(x)={\int }_{{\phi }_1(x)}^{{\phi }_2(x)}f(x,y)\mathrm{d}y,$$

$$I={\int }_a^{b}F(x)\mathrm{d}x\approx \sum _{i=0}^n{\omega }_i^{(1)}F(x_i).$$

为了求出$F(x_i)$,将区间 $[{\phi }_1(x),{\phi }_2(x)]$ 分为m份. 利用单重积分求积公式得:

$$F(x_i)={\int }_{{\phi }_1(x_i)}^{{\phi }_2(x_i)}f(x_i,y)\mathrm{d}y\approx \sum _{j=0}^m{\omega }_j^{(2)}f(x_i,y_j).$$

${\omega }_i^{(1)}\text{ 和 }{\omega }_j^{(2)}$ 分别是 x 方向和 y 方向的求积系数.记 ${\omega }_{ij}={\omega }_i^{(1)}\cdot {\omega }_j^{(2)}$ ,则有二重积分的求积公式:

$$I={\iint }_{\Omega }f(x,y)\mathrm{d}x\mathrm{d}y\approx \sum _{i=0}^n\sum _{j=0}^m{\omega }_{ij}f(x_i,y_j).\text{\hspace{1em}(8)}$$

以上就是二重积分的思路.

因此一般来说 , 计算高维数值积分用复合辛普森公式以及前几节介绍的方法都失效.

6.2 蒙特卡罗方法

一般地 , 我们将随机数序列的模拟方法称为蒙特卡罗方法 .

记 $X_i(i=1,2,\cdot \cdot \cdot ,N)$ 为相互独立的服从区间 [0,1] 上均匀分布的随机变量,则:

$$I_N[f]=\frac{1}{N}\sum _{i=1}^{N}f(X_i).$$

$X_i$ 可通过 N 次独立实验得到 . 由概率论中的大数定理 , 当 N → ∞ 时 , I N [f] 依概率收敛于 f(x)的积分值, 即:

$$I_N[f]\stackrel{\mathrm{P}}{\to }I[f].$$

可以估计f的方差 $Var(f)$ ,取随机变量X的N个样本值 $X_1,X_2,\cdots ,X_N$ ,令:

$${\overline{I}}_N=\frac{1}{N}\sum _{i=1}^{N}f(X_i).$$

then:

$$\mathrm{Var}(f)\approx \frac{1}{N}\sum _{i=1}^N(f(X_i)-{\overline{I}}_N{)}^2$$

这个思路推广到高维积分:

$$I={\int }_0^1\cdots {\int }_0^{1}f(x_1,x_2,\cdots ,x_d)\mathrm{d}{x}_1\mathrm{d}{x}_2\cdots \mathrm{d}{x}_d.$$

仍可取 N 次独立样本值的算术平均:

$$I_N[f]=\frac{1}{N}\sum _{i=1}^{N}f(x_1^{(i)},x_2^{(i)},\cdots ,x_d^{(i)})$$

作为多重积分 I 的无偏估计量 , 其方差仍为$\frac{\mathrm{Var}(f)}{N}$,收敛速度不受积分区域维数 d 的影响.

在实际计算中,$X_i$一般是程序生成的伪随机数.matlab的rand() 可以产生服从均匀分布U(0,1) 的随机数 ( 实际上是伪随机数 ).python的numpy.random.rand() 可以产生服从均匀分布U(0,1) 的随机数.

• 推广到[a,b]的积分范围

那么,如果计算积分:

$$I={\int }_a^{b}f(x)dx$$

对应的蒙特卡洛积分为

$$I^N=\frac{1}{N}\sum _{i=1}^N\frac{f(X_i)}{pdf(X_i)}$$

其中, $pdf(X_i)$ 是概率密度函数, 它描述了随机变量X的概率分布.之前的X_{U(0,1),pdf(X)=1/(1-0),如果X}U(a,b),则pdf(X)=1/(b-a).

• 积分范围推广到[-inf,+inf]的蒙特卡洛方法

积分:

$$F={\int }_{-\infty }^{+\infty }f(x)\mathrm{d}x.$$

可以看成是服从于正态分布随机变量 X 的函数 h(X) 的数学期望, 通过 N 次独立取样的算术平均近似计算:

$$F\approx F_N=\frac{1}{N}\sum _{i=1}^{N}h(x_i),h(x)=\sqrt{2\pi }{\mathrm{e}}^{\frac{x^2}{2}}f(x).$$

可以证明:$E(F_N)=F$

资料链接:

• https://wyman1024.github.io/monte-carlo/
• https://cosx.org/2010/03/monte-carlo-method-to-compute-integration
• https://reichtumqian.github.io/ScientificComputingWiKi/

数值微分

讨论如下问题：已知函数 f(x) 在离散点处的函数值 f(x i )(i = 1,2,··· ,n), 求 f(x)导数.

介绍两种近似方法：

• 基于拉格朗日插值多项式的求导方法
• 基于样条函数的求导方法

7.基于拉格朗日插值多项式的求导方法

已知 $f(x)$ 在离散点处的函数值 $f(x_i)(i=1,2,\cdot \cdot \cdot ,n)$, 可以作出 n 次拉格朗日插值多项式 $L_n(x)$ ,再利用 $L_n(x)$ 的导数来近似代替 $f(x)$ 的导数 , 此方法称为基于拉格朗日插值多项式的数值求导方法

插值余项为:

$$R_n(x)=f(x)-L_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}\Pi (x),\xi \in (a,b).$$

$$R_n^{\prime }(x)=f^{\prime }(x_i)-L_n^{\prime }(x_i)=\frac{f^{(n+1)}(\xi )}{(n+1)!}{\Pi }^{\prime }(x_i)=\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{j=0,j\ne i}^n(x_i-x_j).\text{\hspace{1em}(10)}$$

基于拉格朗日插值的求导方法的缺点是

• 只能求出节点处的导数
• 数据 $f(x_j)$ 有一定误差以及h很小时,数值误差较大,即算法不稳定.
• 高次插值会产生龙格现象,因此实际应用中多采用低次拉格朗日插值型求导公式

以下是基于等距节点的低阶求导计算公式

两点公式

作 f(x) 的线性插值公式:

$$L_1(x)=\frac{x-x_1}{x_0-x_1}f(x_0)+\frac{x-x_0}{x_1-x_0}f(x_1),$$

求导后,令$x=x_0,x_1-x_0=h,$得:

$$L_1^{\prime }(x_0)=\frac{f(x_1)-f(x_0)}{h}.$$

根据公式(10),余项为 $\frac{f^{\prime \prime }(\xi )}{2}(x_0-x_1)$,即:

$$f^{\prime }(x_0)=\frac{f(x_1)-f(x_0)}{h}-\frac{h}{2}{f}^{\prime \prime }(\xi ),\xi \in (x_0,x_1).\text{\hspace{1em}(11-1)}$$

$$f^{\prime }(x_1)=\frac{f(x_1)-f(x_0)}{h}+\frac{h}{2}{f}^{\prime \prime }(\xi ),\xi \in (x_0,x_1).\text{\hspace{1em}(11-2)}$$

略去余项,(11-1),(11-2)称为计算导数的两点公式.分别是向前差商公式和向后差商公式,精度1阶,误差$O(h)$.

三点及五点公式

两个公式共同点在于,在一些点处的微分算不出来,微分个数损失.https://www.math.pku.edu.cn/teachers/lidf/docs/statcomp/html/_statcompbook/approx-diff.html

作 f(x) 的二次插值多项式:

$$L_2(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\frac{(x-x_0)(x-x_1)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2).$$

令 $x=x_0+th,x_i=x_0+ih,$ 上式变成:

$$L_2(x_0+th)=\frac{1}{2}(t-1)(t-2)f(x_0)-t(t-2)f(x_1)+\frac{1}{2}t(t-1)f(x_2),$$

求导,并令t=0,1,2,得:

$$\begin{array}{rl}& f^{\prime }(x_0)=\frac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)]+\frac{h^2}{3}{f}^{(3)}(\xi ),\\ & f^{\prime }(x_1)=\frac{1}{2h}[-f(x_0)+f(x_2)]-\frac{h^2}{6}{f}^{(3)}(\xi ),\\ & f^{\prime }(x_2)=\frac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)]+\frac{h^2}{3}{f}^{(3)}(\xi ).\end{array}\text{\hspace{1em}(12)}$$

以上是带余项的计算导数的三点公式,精度2阶,误差$O(h^2)$.(12-2)为中心差商公式.

同理,五点求导公式为:

$$\begin{array}{rl}& f^{\prime }(x_0)=\frac{1}{12h}(-25f_0+48f_1-36f_2+16f_3-3f_4)+\frac{h^4}{5}{f}^{(5)}(\xi ),\\ & f^{\prime }(x_1)=\frac{1}{12h}(-3f_0-10f_1+18f_2-6f_3+f_4)-\frac{h^4}{20}{f}^{(5)}(\xi ),\\ & f^{\prime }(x_2)=\frac{1}{12h}(f_0-8f_1+8f_3-f_4)-\frac{h^4}{30}{f}^{(5)}(\xi ),\\ & f^{\prime }(x_3)=\frac{1}{12h}(-f_0+6f_1-18f_2+10f_3+3f_4)-\frac{h^4}{20}{f}^{(5)}(\xi ),\\ & f^{\prime }(x_4)=\frac{1}{12h}(3f_0-16f_1+36f_2-48f_3+25f_4)+\frac{h^4}{5}{f}^{(5)}(\xi ).\end{array}\text{\hspace{1em}(13)}$$

五点公式一般比三点公式精确 , 三点公式一般比两点公式精确 ;. 但在实际计算中 , 由于数据 $f_i$ 有误差 , 并不是h 越小计算效果越好

8. 基于样条函数的求导方法

本节我们将用样条函数 $S(x)$ 代替拉格朗日插值多项式 $L_n(x)$作为函数 $f(x)$的近似.若$f(x)\in C^4[a,b],$则有估计:

$$\underset{a⩽x⩽b}{\max }|f^{(k)}(x)-S^{(k)}(x)|⩽C{h}^{4-k},k=0,1,2,3,$$

其中,$h={\max }_{0⩽i⩽n-1}{h}_i.$

此时不但 f(x) 与 S(x) 的函数值很 “ 接近 ”, 它们的导数值也很 “ 接近 ”

设$a=x_0<x_1<\cdots <x_n=b,h_i=x_{i+1}-x_i(i=0,1,\cdots ,n-1).$,且样条函数S(x)在$x_i$处的导数为$S^{\prime }(x_i)=m_i$,则在区间$(x_i,x_{i+1})$上,S(x)为三次多项式,即:

$$\begin{array}{c}{S}_i(x)=\frac{h_i+2(x-x_i)}{h_i^3}(x-x_{i+1}{)}^2{f}_i+\frac{h_i-2(x-x_{i+1})}{h_i^3}(x-x_i{)}^2{f}_{i+1}\\ +\frac{(x-x_i)(x-x_{i+1}{)}^2}{h_i^2}{m}_i+\frac{(x-x_{i+1})(x-x_{i+1}{)}^2}{h_i^2}{m}_{i+1}.\end{array}\text{\hspace{1em}(14)}$$

想要求出$m_i(i=0,1,\cdots ,n)$.需要$S^{\prime \prime }(x)$在$x_i$处连续以及边界条件.

1. 若$S^{\prime }(x_0)=f_0^{\prime }=m_0,S^{\prime }(x_n)=f_n^{\prime }=m_n$已知.m_i满足方程组:

上述方程组 (5.62), (5.65), (5.66) 称为三转角方程组 . 可用追赶法进行求解 . 再将计算结果代入式 (5.61) 中便得 S(x) 以及 S(x) 在任意点 x 处的导数值

由于上述求导数的方法需要求解线性方程组 , 因此也称为隐式格式 . 它具有数值计算稳定的特点 , 可以求出任意点 ( 不一定为节点 ) 的导数值 , 而且精度较高 .

数值实验

# 数值实验t1
from formu_lib import *
def Com3pGausslengendre(f,a:float,b:float,n:int):# 复合三点高斯-勒让德计算积分x=[a+i*(b-a)/n for i in range(n+1)]return sum([Integ1dGuassLegendre(f,x[i],x[i+1],3)for i in range(n)])# 积分函数
pi=lambda x: 4/(1+x**2)# 积分范围
a,b=0,1# 复合simpson计算积分
r1=VarStepInteg1d(pi,a,b,1e-8,method=3)r2=Com3pGausslengendre(pi,a,b,100)
showError(r1,np.pi)
showError(r2,np.pi)
# 数值解: 3.1415926535528365, 精确解: 3.141592653589793, 误差: 1.1763670223853885e-11
# 数值解: 3.1415926535897927, 精确解: 3.141592653589793, 误差: 1.4135798584282297e-16

# t2
from formu_lib import *# 积分函数
f=lambda x: np.sqrt(4-np.sin(x)**2)
a,b=0,np.pi/4r1=VarStepInteg1d(f,a,b,1e-6,method=3)
showError(r1,1.53439197)
# 数值解: 1.5343920054108953, 精确解: 1.53439197, 误差: -2.3078128678621626e-08

求解代码:

from formu_lib import *# 积分函数
f8=lambda x: 1/(1+x**2)
a,b=-5,5result=[]
accurateVal=[]
ns=[]
for n in range(1,51):result.append(Integ1dNewtonCotes(f8,a,b,n))ns.append(n)accurateVal.append(2*np.arctan(5))from matplotlib import pyplot as plt
plt.plot(ns,result,label='Numerical solution')
plt.plot(ns,accurateVal,label='Accurate solution')
plt.xlabel('Newton-Cotes points number')
plt.ylabel(' Integrate Value')
plt.title('Newton-Cotes Integration Value VS Points Number')
plt.show()

从上图可以看出,当随着n的增加, 牛顿-科特斯积分的逐渐不稳定, 但误差也越来越大.不收敛.

code:

# t4
from formu_lib import *
x=[1900,1910,1920,1930,1940,1950,1960,1970,1980,1990]
fs=[76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4]df1=diff1dForward(x,fs)
df2=diff1dBackward(x,fs)
df3=diff1dCentral(x,fs)
plotLines([x,x[:],x[:],x],[fs,df1,df2,df3],["Population","Population growth rate-forward","Population growth rate-backward","Population growth rate-center"])

资料:显式和隐式方法 - 维基百科，自由的百科全书

取步长为h,则: $x_i=a+ih,i=0,1,2,...,n,h=\frac{b-a}{n},x_{i+1}-x_i=h.$ 使用向前差商公式,则有:

$$y^{\prime }(x_i)=\frac{y(x_{i+1})-y(x_i)}{h}$$

代入方程后:

$$y(x_{i+1})-y(x_i)=hf[x_i,y(x_i)],$$

$$y(x_{i+1})=y(x_i)+hf[x_i,y(x_i)],$$

这是显式格式.

如果使用向后差商公式,则有:

$$y^{\prime }(x_{i+1})=\frac{y(x_{i+1})-y(x_i)}{h}$$

代入方程后:

$$y(x_{i+1})-y(x_i)=hf[x_{i+1},y(x_{i+1})],$$

$$y(x_{i+1})-hf[x_{i+1},y(x_{i+1})]=y(x_i)$$

这是隐式格式

• 验证

取$f(x,y)=-y^2,a=0,y(0)=1,b=5>x>a,n取值为100$,则:

# t6
from formu_lib import *
import numpy as np
a,b,n=0,5,100
h=(b-a)/n
x=[a+i*h for i in range(n+1)]f=lambda x: -x*x# result var
ye,yi=np.zeros(n+1),np.zeros(n+1)
ye[0],yi[0]=1,1# explicit method
for i in range(n):ye[i+1]=ye[i]+h*f(ye[i])# implicit method
for i in range(n):yi[i+1]=(-1-np.sqrt(1+4*h*yi[i]))/(2*h)plotLines([x,x],[ye,yi],["y(x)-explictSol","y(x)-implicitSol"])

精确结果:

数值积分微分部分完结