精华内容
下载资源
问答
  • 深度学习求解偏微分方程

    千次阅读 2021-05-25 21:33:21
    深度学习求解偏微分方程1. 稀疏回归解偏微分方程2. 离散连续方程解偏微分方程3. 物理神经网络解偏微分方程(PINN:物理激发的神经网络) 1. 稀疏回归解偏微分方程 论文:《Data-driven discovery of partial ...

    1. 稀疏回归解偏微分方程

    论文:《Data-driven discovery of partial differential equations》
    作者:Samuel H.Rudy, Steven L.Brunton
    具体操作:对同一个变量,用遍历算法逐个筛选出来,所筛选出来的和训练数据对上的即为所求的算子。
    案例:想求解x,y,w之间的关系,就把这三个数据推带入到PDE-FIND算法中,首先找出他们的关系,比如∂w/∂x=C ∂w/∂y 让迭代的训练误差最小,输出的即为稀疏回归的解,然后去拟合参数,拟合出的结果可能是∂w/∂x=0.8 ∂w/∂y
    缺点:想找出训练数据对应的控制方程,就需要先知道原方程的大体样子。通用性不强

    2. 离散连续方程解偏微分方程

    论文:《Learning data0driben discretization for partial differential equations》
    作者:Yohai Bar-Sinai, Stephan Hoyer(Google AI团队)
    具体步骤:先利用传统方法对一个特定的方程进行离散获得高精度的数据,通过训练这些数据来学习如何利用这些离散数据来求近似值导数。
    案例:大自然的景色好比原方程表述的物理信息,用相机去拍照就是用像素点对自然景色进行离散化和逼近,而在方程中,就是用神经网络学习高精度数值解,可以学到微分导数算子,这就类似通过高分辨照片可以学到自然景色一样。
    优点:此方程通过学习离散的数据捕捉到“藏在数据下面”的物理信息,准确度比传统的PDE离散方程高得多,更是省去了传统方法中需要对PDE构造不同的复杂离散格式。

    3. 物理神经网络解偏微分方程(PINN:物理激发的神经网络)

    论文:《Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations》
    作者:M.Raissi, P.Perdikaris
    具体操作:所谓的物理神经网络,其实就是把物理方程作为限制加入神经网络中使训练的结果满足物理规律,即通过把物理方程的迭代前后的差值加到神经网络的损失函数里,让物理方程也参与到了训练过程。
    案例:博格斯方程∂u/∂t+u ∂u/∂x=0.01/π (∂^2 u)/(∂^2 x),u = neural_net(tf.concat([t,x],1), weights, biases),u_t = tf.gradients(u, t)[0],u_x = tf.gradients(u, x)[0],u_xx = tf.gradients(u_x, x)[0]f = u_t + uu_x - (0.01/tf.pi)u_xx。一直训练网络让f和u的损失函数的和迭代不断减小到一定值后,网络就训练好了,就能“练出来”想要找到的函数关系。

    展开全文
  • Python小白的数学建模课-11.偏微分方程数值解法

    千次阅读 多人点赞 2021-08-17 14:10:06
    偏微分方程可以描述各种自然和工程现象, 是构建科学、工程学和其他领域的数学模型主要手段。 本文采用有限差分法求解偏微分方程,通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆...

    偏微分方程可以描述各种自然和工程现象, 是构建科学、工程学和其他领域的数学模型主要手段。

    偏微分方程主要有三类:椭圆方程,抛物方程和双曲方程。

    本文采用有限差分法求解偏微分方程,通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程等常见类型的偏微分方程的数值解法,给出了全部例程和运行结果。

    欢迎关注『Python小白的数学建模课 @ Youcans』系列,每周持续更新。





    1. 偏微分方程基本知识

    微分方程是指含有未知函数及其导数的关系式,偏微分方程是包含未知函数的偏导数(偏微分)的微分方程。

    偏微分方程可以描述各种自然和工程现象, 是构建科学、工程学和其他领域的数学模型主要手段。 科学和工程中的大多数实际问题都归结为偏微分方程的定解问题,如:波传播,流动和扩散,振动,固体力学,电磁学和量子力学,等等。

    偏微分方程主要有三类:椭圆方程,抛物方程和双曲方程。

    • 双曲方程描述变量以一定速度沿某个方向传播,常用于描述振动与波动问题。
    • 椭圆方程描述变量以一定深度沿所有方向传播,常用于描述静电场、引力场等稳态问题。
    • 抛物方程描述变量沿下游传播,常用于描述热传导和扩散等瞬态问题。

    偏微分方程的定解问题通常很难求出解析解,只能通过数值计算方法对偏微分方程的近似求解。常用偏微分方程数值解法有:有限差分方法、有限元方法、有限体方法、共轭梯度法,等等。通常先对问题的求解区域进行网格剖分,然后将定解问题离散为代数方程组,求出在离散网格点上的近似值。

    Python 语言求解偏微分方程的功能是比较弱的,主要有 Fipy, FEniCS 等有限元方法的工具包,另外还有机器学习工具如 Tensorflow 也可以进行偏微分方程的仿真模拟。但是,这些工具都不适合 Python 小白学习和使用,因此本篇采用比较简单的有限差分法对 5种典型的偏微分方程进行编程,通过案例讲解偏微分方程的数值解法。



    2. 案例一:一维线性平流方程

    2.1 一维线性平流方程的数学模型

    平流过程是大气运动中重要的过程。平流方程(Advection equation)描述某一物理量的平流作用而引起局地变化的物理过程,最简单的形式是一维平流方程。

    { ∂ u ∂ t + v ∂ u ∂ x = 0 u ( x , 0 ) = F ( x ) 0 ≤ t ≤ t e x a < x < x b \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x}= 0\\ &u(x,0)=F(x)\\ &0 \leq t \leq t_e\\ &x_a<x<x_b \end{aligned} \end{cases} tu+vxu=0u(x,0)=F(x)0ttexa<x<xb

    式中 u 为某物理量,v 为系统速度,x 为水平方向分量,t 为时间。

    虽然该方程可以求得解析解:
    u ( x , t ) = F ( x − v ∗ t ) ,   0 ≤ t ≤ t e u(x,t)=F(x-v*t), \ 0 \leq t \leq t_e u(x,t)=F(xvt), 0tte

    考虑一维线性平流偏微分方程的数值解法,采用有限差分法求解。简单地, 采用一阶迎风格式的差分方法(First-order Upwind),一阶导数的差分表达式为:
    ( ∂ u ∂ x ) i , j = u i + 1 , j − u i , j Δ x + O ( Δ x ) (\frac{\partial u}{\partial x})_{i,j}=\frac{u_{i+1,j}-u_{i,j}}{\Delta x}+O(\Delta x) (xu)i,j=Δxui+1,jui,j+O(Δx)
    于是得到差分方程:
    u i , j + 1 = u i , j − v ∗ d t / d x ∗ ( u i , j − u i − 1 , j ) u_{i,j+1}=u_{i,j}-v*dt/dx*(u_{i,j}-u_{i-1,j}) ui,j+1=ui,jvdt/dx(ui,jui1,j)
    即可递推求得该平流方程的数值解。


    2.2 偏微分方程编程步骤

    以该题为例一类有限差分法求解一维平流问题偏微分方程的步骤:

    1. 导入numpy、matplotlib 包;
    2. 定义初始条件函数 U(x,0);
    3. 输入模型参数 v, p,定义求解的时间域 (tStart, tEnd) 和空间域 (xMin, xMax),设置差分步长 dt, nNodes;
    4. 初始化;
    5. 递推求解差分方程在区间 [xa,xb] 的数值解,获得网格节点的处的 u(t) 值。

    在例程中,设初值条件为 F ( x ) = s i n ( 2 ∗ ( x − p ) 2 ) F(x) = sin(2 * (x-p)^2) F(x)=sin(2(xp)2),取 v = 1.0 , p = 0.25 v= 1.0, p=0.25 v=1.0,p=0.25,空间域 x ∈ ( 0 , π ) x\in(0,\pi) x(0,π)


    2.3 Python 例程:偏微分方程(一维平流方程)

    # mathmodel13_v1.py
    # Demo10 of mathematical modeling algorithm
    # Solving partial differential equations
    # 偏微分方程数值解法
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    # 1. 一维平流方程 (advection equation)
    # U_t + v*U_x = 0
    
    # 初始条件函数 U(x,0)
    def funcUx0(x, p): 
        u0 = np.sin(2 * (x-p)**2)
        return u0
    
    # 输入参数
    v1 = 1.0  # 平流方程参数,系统速度
    p = 0.25  # 初始条件函数 u(x,0) 中的参数
    
    tc = 0  # 开始时间
    te = 1.0  # 终止时间: (0, te)
    xa = 0.0  # 空间范围: (xa, xb)
    xb = np.pi
    dt = 0.02  # 时间差分步长
    nNodes = 100  # 空间网格数
    
    # 初始化
    nsteps = round(te/dt)
    dx = (xb - xa) / nNodes
    x = np.arange(xa-dx, xb+2*dx, dx)
    ux0 = funcUx0(x, p)
    
    u = ux0.copy()  # u(j)
    ujp = ux0.copy()  # u(j+1)
    
    # 时域差分
    for i in range(nsteps):
        plt.clf()  # 清除当前 figure 的所有axes, 但是保留当前窗口
    
        # 计算 u(j+1)
        for j in range(nNodes + 2):
            ujp[j] = u[j] - (v1 * dt/dx) * (u[j] - u[j-1])
    
        # 更新边界条件
        u = ujp.copy()
        u[0] = u[nNodes + 1]
        u[nNodes+2] = u[1]
    
        # 绘图
        plt.plot(x, u, 'b-', label="v1= 1.0")
        plt.axis((xa-0.1, xb + 0.1, -1.1, 1.1))
        plt.xlabel("x")
        plt.ylabel("U(x)")
        plt.legend(loc=(0.05,0.05))
        plt.title("Advection equation with finite difference method, t = %1.f" % (tc + dt))
        plt.text(0.05,0.9,"youcans-xupt",color='gainsboro')
        plt.pause(0.001)
        tc += dt
    
    plt.show()
    

    2.4 Python 例程运行结果

    在这里插入图片描述



    3. 案例二:一维热传导方程

    3.1 一维热传导方程的数学模型

    热传导方程描述一个区域内的温度随时间的变化,是典型的抛物型偏微分方程,也称为扩散方程。

    一维热传导方程考虑各向同性的均匀细杆,在垂直于 x 轴的截面上的温度相同,细杆的圆周与周围环境无热交换,杆内无热源,则温度 u = u ( t , x ) u=u(t,x) u=u(t,x) 是时间变量 t 和水平方向空间变量 x 的函数。

    { ∂ u ∂ t = d i v ( U u ) = λ ∂ 2 u ∂ x 2 u ( x , 0 ) = u 0 ( x ) u ( x a ) = u a ( t ) ,   u ( x b , t ) = u b ( t ) 0 ≤ t ≤ t e ,   x a < x < x b \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} = div(U_u) = \lambda \frac{\partial ^2u}{\partial x^2}\\ &u(x,0)=u_0(x)\\ &u(x_a)=u_a(t),\ u(x_b,t)=u_b(t)\\ &0\leq t \leq t_e, \ x_a < x < x_b \end{aligned} \end{cases} tu=div(Uu)=λx22uu(x,0)=u0(x)u(xa)=ua(t), u(xb,t)=ub(t)0tte, xa<x<xb

    式中 λ \lambda λ 为热扩散率,取决于材料本身的热传导率、密度和热容。

    考虑一维热传导偏微分方程的数值解法,采用有限差分法求解。简单地, 采用迎风法的三点差分格式,二阶导数的差分表达式为:
    ( ∂ 2 u ∂ x 2 ) i , j = u i + 1 , j − 2 u i , j + u i − 1 , j ( Δ x ) 2 + O ( Δ x ) 2 (\frac{\partial ^2u}{\partial x^2})_{i,j} = \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}+O(\Delta x)^2 (x22u)i,j=(Δx)2ui+1,j2ui,j+ui1,j+O(Δx)2
    于是得到差分方程:
    u i , j + 1 = u i , j + λ d t / d x 2 ∗ ( u i + 1 , j − 2 u i , j + u i − 1 , j ) u_{i,j+1} = u_{i,j} + \lambda dt/dx^2*(u_{i+1,j}-2u_{i,j}+u_{i-1,j}) ui,j+1=ui,j+λdt/dx2(ui+1,j2ui,j+ui1,j)
    即可递推求得一维热传导方程的数值解。


    3.2 偏微分方程编程步骤

    以该题为例一类有限差分法求解一维平流问题偏微分方程的步骤:

    1. 导入numpy、matplotlib 包;
    2. 输入模型参数 L, lam,定义求解的时间域 (0, te) 和空间域 (0, L),设置差分步长 dt, dx;
    3. 初始化;
    4. 计算每一时点的边界条件 U[0,j],U[L,j]、每一位置的初始条件U[i,0];
    5. 递推求解差分方程在空间域 [0, L] 的数值解,获得网格节点的处的 U(x,t) 。

    在例程中,设初值条件为 u ( x , t = 0 ) = 0 u(x,t=0) = 0 u(x,t=0)=0,边界条件为 u ( x a , t ) = F ( t ) ,   u ( x b , t ) = 0 u(x_a,t)=F(t), \ u(x_b,t)=0 u(xa,t)=F(t), u(xb,t)=0,在时间域 t ∈ ( 0 , 2.0 ) t\in(0,2.0) t(0,2.0)、空间域 x ∈ ( 0 , 1.0 ) x\in(0,1.0) x(0,1.0) 求数值解即温度分布。
    (1)例程中的初始条件 U[i,0] 为常数,如果初始条件是 x 的函数 u(x,0),将该函数在初始条件语句赋值即可(参加例程中注释的语句)。
    (2)例程中的边界条件,一端是 t 的函数 u(0,t),另一端是常数 u(L,t) =0,这些条件也可以根据具体问题设置为相应的常数或函数。


    3.3 Python 例程:偏微分方程(一维热传导方程)

    # mathmodel13_v1.py
    # Demo10 of mathematical modeling algorithm
    # Solving partial differential equations
    # 偏微分方程数值解法
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    # 2. 一维热传导方程(抛物型偏微分方程)
    # pu/pt = l*p2u/px2
    
    # 模型参数
    L = 1.0  # 细杆长度
    lam = 1.0  # 热扩散率
    tc = 0  # 开始时间
    te = 10.0  # 终止时间: (0, te)
    
    # 初始化
    dx = 0.05  # 空间步长
    dt = 0.001  # 时间步长
    nNodes = round(L/dx)  # 空间网格数
    nSteps = round(te/dt)  # 时序网格数
    K = lam * dt/(dx**2)  # lambda * dt/dx^2
    U = np.zeros([nNodes+1, nSteps+1])  # 建立二维数组
    
    # 边界条件
    for j in np.arange(0, nSteps+1):  # 时间序列
        U[0,j] = 7.5 + (nSteps-j)/2000 * np.sin(j/1000)/(1+np.exp(-j))
        U[nNodes,j] = 0.0  # 每一时点的边界条件
    
    # 初始条件
    for i in np.arange(0, nNodes):  # 空间序列
        # U[i,0]= 0.2*i*h*(L-i*h)  # 初始条件是 x 的函数
        U[i,0]= 0  # 每一位置的初始条件
    
    # 时域差分法求解
    for j in np.arange(0, nSteps):  # 时间步长
        for i in np.arange(1, nNodes):  # 空间步长
            U[i,j+1] = K*U[i+1,j] + (1-2*K)*U[i,j] + K*U[i-1,j]
    
    # 绘图
    xZone = np.arange(0, (nNodes+1)*dx, dx)  # 建立空间网格
    tZone = np.arange(0, (nSteps+1)*dt, dt)  # 建立空间网格
    fig = plt.figure(figsize=(10, 6))
    rect1 = [0.05, 0.2, 0.4, 0.65]  # [左, 下, 宽, 高], 0.0~1.0
    ax1 = plt.axes(rect1)
    for k in range(0,nSteps+1,round(nSteps/5)):
        ax1.plot(xZone, U[:,k], label=r"x={}".format(k/nSteps))
    ax1.set_ylabel('u(x,t)')
    ax1.set_xlabel('x')
    ax1.set_xlim(0,L)
    ax1.set_ylim(-1,12)
    ax1.set_title("Temperature distribution along t-axis")
    ax1.legend(loc='upper right')
    rect2 = [0.55, 0.2, 0.4, 0.65]  # [左, 下, 宽, 高], 0.0~1.0
    ax2 = plt.axes(rect2)
    for k in range(0,nNodes+1,round(nNodes/5)):  # U[nNodes,k] = 0.0
        ax2.plot(tZone, U[k,:], label=r"t={}".format(k/nNodes))
    ax2.set_ylabel('u(x,t)')
    ax2.set_xlabel('t')
    ax2.set_xlim(0,te)
    ax2.set_ylim(-1,12)
    ax2.set_title("Temperature distribution along x-axis")
    ax2.legend(loc='upper right')
    plt.show()
    

    3.4 Python 例程运行结果

    在这里插入图片描述



    4. 案例三:二维双曲型方程

    4.1 二维波动方程的数学模型

    波动方程(wave equation)是典型的双曲形偏微分方程,广泛应用于声学,电磁学,和流体力学等领域,描述自然界中的各种的波动现象,包括横波和纵波,例如声波、光波和水波。

    考虑如下二维波动方程的初边值问题:
    { ∂ 2 u ∂ t 2 = c 2 ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) ∂ ∂ t u ( 0 , x , y ) = 0 u ( x , y , 0 ) = u 0 ( x , y ) u ( 0 , y , t ) = u a ( t ) ,   u ( 1 , y , t ) = u b ( t ) u ( x , 0 , t ) = u c ( t ) ,   u ( x , 1 , t ) = u d ( t ) 0 ≤ t ≤ t e ,   0 < x < 1 ,   0 < y < 1 \begin{cases} \begin{aligned} &\frac{\partial^2 u}{\partial t^2} = c^2 (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2})\\ &\frac{\partial }{\partial t} u(0,x,y)= 0\\ &u(x,y,0)=u_0(x,y)\\ &u(0,y,t)=u_a(t),\ u(1,y,t)=u_b(t)\\ &u(x,0,t)=u_c(t),\ u(x,1,t)=u_d(t)\\ &0\leq t \leq t_e, \ 0 < x < 1, \ 0< y < 1 \end{aligned} \end{cases} t22u=c2(x22u+y22u)tu(0,x,y)=0u(x,y,0)=u0(x,y)u(0,y,t)=ua(t), u(1,y,t)=ub(t)u(x,0,t)=uc(t), u(x,1,t)=ud(t)0tte, 0<x<1, 0<y<1

    式中:u 是振幅;c 为波的传播速率,c 可以是固定常数,或位置的函数 c(x,y),也可以是振幅的函数 c(u)。

    考虑二维波动偏微分方程的数值解法,采用有限差分法求解。简单地, 采用迎风法的三点差分格式, 将上述的偏微分方程离散为差分方程 :
    在这里插入图片描述
    化简后得到:
    在这里插入图片描述

    即可递推求得二维波动方程的数值解。

    为了保证算法的收敛性,迎风法的三点差分格式要求步长比小于 1:

    r = 4 c 2 Δ t 2 Δ x 2 + Δ y 2 ≤ 1 r = \frac{4c^2\Delta t^2}{\Delta x^2 + \Delta y^2} \leq 1 r=Δx2+Δy24c2Δt21


    4.2 偏微分方程编程步骤

    以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤:

    1. 导入numpy、matplotlib 包;
    2. 输入模型参数 c,定义求解的时间域 (0, te) 和空间域 (0,1);
    3. 初始化,设置差分步长 dt, dx, dy,检验步长比以保证算法的稳定性;
    4. 计算初值条件U[0]、U[1];
    5. 递推求解差分方程在空间域 [0, 1]*[0, 1] 的数值解,获得网格节点的处的 U(t,x,y) 。

    在例程中,取初始条件为 u ( x , y , 0 ) = s i n ( 6 π x ) + c o s ( 4 π y ) u(x,y,0)=sin(6\pi x)+cos(4\pi y) u(x,y,0)=sin(6πx)+cos(4πy),边界条件为 u ( 0 , y , t ) = u ( 1 , y , t ) = 0 u(0,y,t)=u(1,y,t)=0 u(0,y,t)=u(1,y,t)=0 u ( x , 0 , t ) = u ( x , 1 , t ) = 0 u(x,0,t)=u(x,1,t)=0 u(x,0,t)=u(x,1,t)=0,在时间域 t ∈ ( 0 , 1.0 ) t\in(0,1.0) t(0,1.0)、空间域 x ∈ ( 0 , 1.0 ) ,   y ∈ ( 0 , 1.0 ) x\in(0,1.0),\ y\in(0,1.0) x(0,1.0), y(0,1.0) 求数值解即振动状态。


    4.3 Python 例程

    # mathmodel13_v1.py
    # Demo10 of mathematical modeling algorithm
    # Solving partial differential equations
    # 偏微分方程数值解法
    
    # 4. 二维波动方程(双曲型二阶偏微分方程)
    # p2u/pt2 = c^2*(p2u/px2+p2u/py2)
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    # 模型参数
    c = 1.0  # 波的传播速率
    tc, te = 0.0, 1.0  # 时间范围,0<t<te
    xa, xb = 0.0, 1.0  # 空间范围,xa<x<xb
    ya, yb = 0.0, 1.0  # 空间范围,ya<y<yb
    
    # 初始化
    c2 = c*c  # 方程参数
    dt = 0.01  # 时间步长
    dx = dy = 0.02  # 空间步长
    tNodes = round(te/dt)  # t轴 时序网格数
    xNodes = round((xb-xa)/dx)  # x轴 空间网格数
    yNodes = round((yb-ya)/dy)  # y轴 空间网格数
    tZone = np.arange(0, (tNodes+1)*dt, dt)  # 建立空间网格
    xZone = np.arange(0, (xNodes+1)*dx, dx)  # 建立空间网格
    yZone = np.arange(0, (yNodes+1)*dy, dy)  # 建立空间网格
    xx, yy = np.meshgrid(xZone, yZone)  # 生成网格点的坐标 xx,yy (二维数组)
    
    # 步长比检验(r>1 则算法不稳定)
    r = 4 * c2 * dt*dt / (dx*dx+dy*dy)
    print("dt = {:.2f}, dx = {:.2f}, dy = {:.2f}, r = {:.2f}".format(dt,dx,dy,r))
    assert r < 1.0, "Error: r>1, unstable step ratio of dt2/(dx2+dy2) ."
    rx = c*c * dt**2/dx**2
    ry = c*c * dt**2/dy**2
    
    # 绘图
    fig = plt.figure(figsize=(8,6))
    ax1 = fig.add_subplot(111, projection='3d')
    # ax2 = fig.add_subplot(122, projection='3d')
    
    # 计算初始值
    U = np.zeros([tNodes+1, xNodes+1, yNodes+1])  # 建立三维数组
    U[0] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[0,:,:]
    U[1] = np.sin(6*np.pi*xx)+np.cos(4*np.pi*yy)  # U[1,:,:]
    surf = ax1.plot_surface(xx, yy, U[0,:,:], rstride=2, cstride=2, cmap=plt.cm.coolwarm)
    # wframe = ax2.plot_wireframe(xx, yy, U[0], rstride=2, cstride=2, linewidth=1)
    
    # 有限差分法求解
    for k in range(2,tNodes+1):
        if surf:
            ax1.collections.remove(surf)  # 更新三维动画窗口
    
        for i in range(1,xNodes):
            for j in range(1,yNodes):
                U[k,i,j] = rx*(U[k-1,i-1,j]+U[k-1,i+1,j]) + ry*(U[k-1,i,j-1]+U[k-1,i,j+1])\
                         + 2*(1-rx-ry)*U[k-1,i,j] -U[k-2,i,j]
    
        surf = ax1.plot_surface(xx, yy, U[k,:,:], rstride=2, cstride=2, cmap='rainbow')
        # wframe = ax2.plot_wireframe(xx, yy, U[k,:,:], rstride=2, cstride=2, linewidth=1, cmap='rainbow')
        ax1.set_xlim3d(0, 1.0)
        ax1.set_ylim3d(0, 1.0)
        ax1.set_zlim3d(-2, 2)
        ax1.set_title("2D wave equationt (t= %.2f)" % (k*dt))
        ax1.set_xlabel("x-youcans")
        ax1.set_ylabel("y-XUPT")
        plt.pause(0.01)
    
    plt.show()
    

    4.4 Python 例程运行结果

    在这里插入图片描述



    5. 案例四:二维抛物型方程

    5.1 二维热传导方程的数学模型

    热传导方程(heat equation)是典型的抛物形偏微分方程,也成为扩散方程。广泛应用于声学,电磁学,和流体力学等领域,描述自然界中的各种的波动现象,包括横波和纵波,例如声波、光波和水波。

    考虑如下二维热传导方程的初边值问题:
    { ∂ u ∂ t = λ ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) + q v u ( x , y , 0 ) = u 0 ( x , y ) u ( 0 , y , t ) = u a ( t ) ,   u ( 1 , y , t ) = u b ( t ) u ( x , 0 , t ) = u c ( t ) ,   u ( x , 1 , t ) = u d ( t ) 0 ≤ t ≤ t e ,   0 < x < 1 ,   0 < y < 1 \begin{cases} \begin{aligned} &\frac{\partial u}{\partial t} = \lambda (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) + q_v\\ &u(x,y,0)=u_0(x,y)\\ &u(0,y,t)=u_a(t),\ u(1,y,t)=u_b(t)\\ &u(x,0,t)=u_c(t),\ u(x,1,t)=u_d(t)\\ &0\leq t \leq t_e, \ 0 < x < 1, \ 0< y < 1 \end{aligned} \end{cases} tu=λ(x22u+y22u)+qvu(x,y,0)=u0(x,y)u(0,y,t)=ua(t), u(1,y,t)=ub(t)u(x,0,t)=uc(t), u(x,1,t)=ud(t)0tte, 0<x<1, 0<y<1

    式中 λ \lambda λ 为热扩散率,取决于材料本身的热传导率、密度和热容; q v q_v qv 是热源强度,可以是恒定值,也可以是时间、空间的函数。

    考虑二维抛物型偏微分方程的数值解法,采用有限差分法求解。将上述的偏微分方程离散为差分方程 :

    在这里插入图片描述

    化简后得到:
    在这里插入图片描述

    即可递推求得二维波动方程的数值解:
    { U k + 1 = U k + r x ∗ A ∗ U k + r y ∗ B ∗ U k + q v ∗ d t ) A = [ − 2 1 ⋯ 0 0 1 − 2 1 ⋯ 0 0 1 − 2 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 − 2 ] ( N x + 1 , N x + 1 ) B = [ − 2 1 ⋯ 0 0 1 − 2 1 ⋯ 0 0 1 − 2 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ 1 − 2 ] ( N y + 1 , N y + 1 ) \begin{cases} \begin{aligned} &U_{k+1} = U_{k} + r_x*A*U_k + r_y*B*U_k + q_v*dt)\\ &A= \begin{bmatrix} -2 & 1 & \cdots & 0 & 0\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots& \vdots\\ 0 & 0 & \cdots & 1 & -2\\ \end{bmatrix} _{(Nx+1,Nx+1)} B= \begin{bmatrix} -2 & 1 & \cdots & 0 & 0\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots& \vdots\\ 0 & 0 & \cdots & 1 & -2\\ \end{bmatrix} _{(Ny+1,Ny+1)} \end{aligned} \end{cases} Uk+1=Uk+rxAUk+ryBUk+qvdt)A=2100121012010002(Nx+1,Nx+1)B=2100121012010002(Ny+1,Ny+1)


    5.2 偏微分方程编程步骤

    以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤:

    1. 导入numpy、matplotlib 包;
    2. 输入热传导参数、热源参数等模型参数,定义求解的时间域 (0, te) 和空间域;
    3. 初始化,设置差分步长 dt, dx, dy,计算三对角系数矩阵 A、B,差分系数 rx, ry, ft;
    4. 计算初始条件 U 0 U_0 U0
    5. 递推求解差分方程在空间域的数值解,更新边界条件,获得网格节点的处的 U ( x , y ) U(x,y) U(x,y)
    6. 绘制等温云图。

    例程求解一个薄板受热的温度分布问题,其初始条件为 t I n i = 25 tIni =25 tIni=25,边界条件为 t B o u n d = 25 tBound = 25 tBound=25,热源为 Q v Qv Qv,在空间域 x ∈ ( 0 , 16 ) ,   y ∈ ( 0 , 12 ) x\in(0,16),\ y\in(0,12) x(0,16), y(0,12) ,时间域 t ∈ ( 0 , 5 ) t\in(0,5) t(0,5) 求数值解即温度分布。

    对于外加热源,例程中给出了三种情况:(1)恒定热源,(2)热源功率(或开关)随时间变化,(3)热源位置随时间变化(从 (x0,y0) 以速度 (xv,yv) 移动,以模拟焊接加热的情况)。


    5.3 Python 例程

    # mathmodel13_v1.py
    # Demo10 of mathematical modeling algorithm
    # Solving partial differential equations
    # 偏微分方程数值解法
    
    # 5. 二维热传导方程(抛物型二阶偏微分方程)
    # pu/p2 = c*(p2u/px2+p2u/py2)
    
    import numpy as np
    import matplotlib.pyplot as plt
    
    def showcontourf(zMat, xyRange, tNow):  # 绘制等温云图
        x = np.linspace(xyRange[0], xyRange[1], zMat.shape[1])
        y = np.linspace(xyRange[2], xyRange[3], zMat.shape[0])
        xx,yy = np.meshgrid(x,y)
        zMax = np.max(zMat)
        yMax, xMax = np.where(zMat==zMax)[0][0], np.where(zMat==zMax)[1][0]
        levels = np.arange(0,100,1)
        showText = "time = {:.1f} s\nhotpoint = {:.1f} C".format(tNow, zMax)
    
        plt.clf()  # 清除当前图形及其所有轴,但保持窗口打开
        plt.plot(x[xMax],y[yMax],'ro')  # 绘制最高温度点
        plt.contourf(xx, yy, zMat, 100, cmap=plt.cm.get_cmap('jet'), origin='lower', levels = levels)
        plt.annotate(showText, xy=(x[xMax],y[yMax]), xytext=(x[xMax],y[yMax]),fontsize=10)
        plt.colorbar()
        plt.xlabel('Xupt')
        plt.ylabel('Youcans')
        plt.title('Temperature distribution of the plate')
        plt.draw()
    
    # 模型参数
    uIni = 25  # 初始温度值
    uBound = 25.0  # 边界条件
    c = 1.0  # 热传导参数
    qv = 50.0  # 热源功率
    x0, y0 = 0.0, 3.0  # 热源初始位置
    vx, vy = 2.0, 1.0  # 热源移动速度
    # 求解范围
    tc, te = 0.0, 5.0  # 时间范围,0<t<te
    xa, xb = 0.0, 16.0  # 空间范围,xa<x<xb
    ya, yb = 0.0, 12.0  # 空间范围,ya<y<yb
    
    # 初始化
    dt = 0.002  # 时间步长
    dx = dy = 0.1  # 空间步长
    tNodes = round(te/dt)  # t轴 时序网格数
    xNodes = round((xb-xa)/dx)  # x轴 空间网格数
    yNodes = round((yb-ya)/dy)  # y轴 空间网格数
    xyRange = np.array([xa, xb, ya, yb])
    xZone = np.linspace(xa, xb, xNodes+1)  # 建立空间网格
    yZone = np.linspace(ya, yb, yNodes+1)  # 建立空间网格
    xx,yy = np.meshgrid(xZone, yZone)  # 生成网格点的坐标 xx,yy (二维数组)
    
    # 计算 差分系数矩阵 A、B (三对角对称矩阵),差分系数 rx,ry,ft
    A = (-2) * np.eye(xNodes+1, k=0) + (1) * np.eye(xNodes+1, k=-1) + (1) * np.eye(xNodes+1, k=1)
    B = (-2) * np.eye(yNodes+1, k=0) + (1) * np.eye(yNodes+1, k=-1) + (1) * np.eye(yNodes+1, k=1)
    rx, ry, ft = c*dt/(dx*dx), c*dt/(dy*dy), qv*dt
    
    # 计算 初始值
    U = uIni * np.ones((yNodes+1, xNodes+1))  # 初始温度 u0
    
    # 前向欧拉法一阶差分求解
    for k in range(tNodes+1):
        t = k * dt  # 当前时间
    
        # 热源条件
        # (1) 恒定热源:Qv(x0,y0,t) = qv, 功率、位置 恒定
        # Qv = qv
        # (2) 热源功率随时间变化 Qv(x0,y0,t)=f(t)
        # Qv = qv*np.sin(t*np.pi) if t<2.0 else qv
        # (3) 热源位置随时间变化 Qv(x,y,t)=f(x(t),y(t))
        xt, yt = x0+vx*t, y0+vy*t  # 热源位置变化
        Qv = qv * np.exp(-((xx-xt)**2+(yy-yt)**2))  # 热源方程
    
        # 边界条件
        U[:,0] = U[:,-1] = uBound
        U[0,:] = U[-1,:] = uBound
    
        # 差分求解
        U = U + rx * np.dot(U,A) + ry * np.dot(B,U) + Qv*dt
    
        if k % 100 == 0:
            showcontourf(U, xyRange, k*dt)  # 绘制等温云图
            print('t={:.2f}s\tTmax={:.1f}  Tmin={:.1f}'.format(t, np.max(U), np.min(U)))
    

    5.4 Python 例程运行结果

    在这里插入图片描述



    6. 案例五:二维椭圆型方程

    6.1 二维椭圆方程的数学模型

    椭圆型偏微分方程是一类重要的偏微分方程,主要用来描述物理的平衡稳定状态,如定常状态下的电磁场、引力场和反应扩散现象等,广泛应用于流体力学、弹性力学、电磁学、几何学和变分法中。
    考虑如下二维泊松方程:
    ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = f ( x , y ) \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)\\ x22u+y22u=f(x,y)

    上式在 f ( x , y ) = 0 f(x,y)=0 f(x,y)=0 时就是拉普拉斯方程(Laplace equation)。

    考虑二维椭圆型偏微分方程的数值解法,采用有限差分法求解。简单地, 采用五点差分格式表示二阶导数的差分表达式,将上述的偏微分方程离散为差分方程 :
    ( u i − 1 , j − 2 u i , j + u i + 1 , j ) / Δ h 2 + ( u i , j − 1 − 2 u i , j + u i , j + 1 ) / Δ h 2 = f i , j (u_{i-1,j}-2u_{i,j}+u_{i+1,j})/\Delta h^2 + (u_{i,j-1}-2u_{i,j}+u_{i,j+1})/\Delta h^2 = f_{i,j} (ui1,j2ui,j+ui+1,j)/Δh2+(ui,j12ui,j+ui,j+1)/Δh2=fi,j
    椭圆型偏微分描述不随时间变化的均衡状态,没有初始条件,因此不能沿时间步长递推求解。对上式的差分方程,可以通过矩阵求逆方法求解,但当 h 较小时网格很多,矩阵求逆的内存占用和计算量极大。于是,可以使用迭代松弛法递推求得二维椭圆方程的数值解:
    u i , j k + 1 = ( 1 − w ) ∗ u i , j k + w / 4 ∗ ( u i , j + 1 k + u i , j − 1 k + u i − 1 , j k + u i + 1 , j k − h 2 ∗ f i , j ) u_{i,j}^{k+1} = (1-w)*u_{i,j}^k + w/4*(u_{i,j+1}^k + u_{i,j-1}^k + u_{i-1,j}^k + u_{i+1,j}^k -h^2* f_{i,j}) ui,jk+1=(1w)ui,jk+w/4(ui,j+1k+ui,j1k+ui1,jk+ui+1,jkh2fi,j)


    6.2 偏微分方程编程步骤

    以该题为例一类有限差分法求解二维波动问题偏微分方程的步骤:

    1. 导入numpy、matplotlib 包;
    2. 输入模型参数,定义空间域 (0,1);
    3. 初始化,设置差分步长 h、松弛因子 w;
    4. 计算边值条件 u[0,:], u[1,:], u[:,0], u[:,1];
    5. 迭代松弛法递推求解差分方程在空间域 [0, 1]*[0, 1] 的数值解 。

    在例程中,取边界条件为 u ( x , 0 ) = u ( x , 1 ) = 0 ,    u ( 0 , y ) = u ( 1 , y ) = 1 u(x,0)=u(x,1)=0, \; u(0,y)=u(1,y)=1 u(x,0)=u(x,1)=0,u(0,y)=u(1,y)=1。为了让图形更有趣,也可以参考例程中选择不同的边界条件。


    6.3 Python 例程

    # mathmodel13_v1.py
    # Demo10 of mathematical modeling algorithm
    # Solving partial differential equations
    # 偏微分方程数值解法
    
    # 7. 二维椭圆方程(Laplace equation)
    # u_xx+u_yy=s
    
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    # 求解范围
    xa, xb = 0.0, 1.0  # 空间范围,xa<x<xb
    ya, yb = 0.0, 1.0  # 空间范围,ya<y<yb
    
    # 初始化
    h = 0.01  # 空间步长, dx = dy = 0.01
    w = 0.5  # 松弛因子
    nodes = round((xb-xa)/h)  # x轴 空间网格数
    
    # 边值条件
    u = np.zeros((nodes+1, nodes+1))
    # u[:, 0] = 0
    # u[:, -1] = 0  # -1 表示数组的最后一个值
    # u[0, :] = 1
    # u[-1, :] = 1  # -1 表示数组的最后一个值
    for i in range(nodes+1):
        u[i, 0] = 1.0 + np.sin(0.5*(i-50)/np.pi)
        u[i, -1] = -1.0 + 0.5*np.sin((i-50)/np.pi)
        u[0, i] = -1.0 + 0.5*np.sin((i-50)/np.pi)
        u[-1, i] = 1.0 + np.sin(0.5*(50-i)/np.pi)
    
    # 迭代松弛法求解
    for iter in range(100):
        for i in range(1, nodes):
            for j in range(1, nodes):
                u[i, j] = w/4 * (u[i-1, j] + u[i+1, j] + u[i, j-1] + u[i, j+1]) + (1-w) * u[i, j]
    
    # 绘图
    x = np.linspace(0, 1, nodes+1)
    y = np.linspace(0, 1, nodes+1)
    xx, yy = np.meshgrid(x, y)
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(xx, yy, u, cmap=plt.get_cmap('rainbow'))
    fig.colorbar(surf, shrink=0.5)
    ax.set_xlim3d(0, 1.0)
    ax.set_ylim3d(0, 1.0)
    ax.set_zlim3d(-2, 2.5)
    ax.set_title("2D elliptic partial differential equation")
    ax.set_xlabel("Xupt")
    ax.set_ylabel("Youcans")
    plt.show()
    

    6.4 Python 例程运行结果

    在这里插入图片描述



    7. 小结

    1. 偏微分方程是数学中的重要学科分支,其数值解法的研究和方法可谓博大精深,远非本文所能涵盖。
    2. 本文针对 Python小白,采用有限差分法求解偏微分方程,通过案例介绍了一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程的数值解法,给出了例程和运行结果,供 Python 小白参考,以便进行数模学习。
    3. 需要特别说明的是:一阶差分格式的精度相对较低,往往需要较多网格计算。随着精度的不断提高,可以推导出各种差分表达式。高阶精度的差分公式可以给出质量更高的解,但需要使用更多的网格点进行计算,程序更复杂,计算量很大。类似地,显式方法的建立和编程简单,但要求步长较小才能保持稳定性,隐式方法的建立和编程更复杂,运算量大,但稳定性好。
    4. 需要特别说明的是:关于偏微分方程数值解法的稳定性、收敛性和误差分析,是非常专业的问题。例如步长选择不恰当可能导致算法不稳定,如 4.2 中应满足步长比 r>1。除非找到专业课程教材或范文中有相关内容可以参考套用,否则不建议小白自己摸索,这些问题不是调整参数试试就能试出来的。

    【本节完】


    版权声明:

    欢迎关注『Python小白的数学建模课 @ Youcans』 原创作品

    原创作品,转载必须标注原文链接:https://blog.csdn.net/youcans/article/details/119755450

    Copyright 2021 Youcans, XUPT

    Crated:2021-08-16



    欢迎关注 『Python小白的数学建模课 @ Youcans』 系列,持续更新
    Python小白的数学建模课-01.新手必读
    Python小白的数学建模课-02.数据导入
    Python小白的数学建模课-03.线性规划
    Python小白的数学建模课-04.整数规划
    Python小白的数学建模课-05.0-1规划
    Python小白的数学建模课-06.固定费用问题
    Python小白的数学建模课-07.选址问题
    Python小白的数学建模课-09.微分方程模型
    Python小白的数学建模课-10.微分方程边值问题
    Python小白的数学建模课-11.偏微分方程数值解法
    Python小白的数学建模课-12.非线性规划
    Python小白的数学建模课-15.图论的基本概念
    Python小白的数学建模课-16.最短路径算法
    Python小白的数学建模课-17.条件最短路径算法
    Python小白的数学建模课-18.最小生成树问题
    Python小白的数学建模课-19.网络流优化问题
    Python小白的数学建模课-20.网络流优化案例
    Python小白的数学建模课-21.关键路径法
    Python小白的数学建模课-22.插值方法
    Python小白的数学建模课-23.数据拟合全集
    Python小白的数学建模课-A1.国赛赛题类型分析
    Python小白的数学建模课-A2.2021年数维杯C题探讨
    Python小白的数学建模课-A3.12个新冠疫情数模竞赛赛题及短评
    Python小白的数学建模课-B2. 新冠疫情 SI模型
    Python小白的数学建模课-B3. 新冠疫情 SIS模型
    Python小白的数学建模课-B4. 新冠疫情 SIR模型
    Python小白的数学建模课-B5. 新冠疫情 SEIR模型
    Python小白的数学建模课-B6. 新冠疫情 SEIR改进模型
    Python数模笔记-PuLP库
    Python数模笔记-StatsModels统计回归
    Python数模笔记-Sklearn
    Python数模笔记-NetworkX
    Python数模笔记-模拟退火算法

    展开全文
  • 前言 在完成了常微分的数值解之后,我开始如法炮制的来... 今天选择的偏微分方程为Burgers方程,在Dirichlet条件下求数值解,具体形式为: 训练流程 1.定义网络 这里定义网络的方法还是和第二篇是一样的,有...

    前言

            在完成了常微分的数值解之后,我开始如法炮制的来解偏微分,我觉得解法上是一样的,都直接使用autograd就可以了,所以理论是难度并不大(虽然实际上我是花的时间最长的),只不过需要注意的细节比较多,所以模型卡了很久都没有效果,让我都怀疑人生了,今天终于找到了问题!所以把整个过程介绍一下

            今天选择的偏微分方程为Burgers方程,在Dirichlet条件下求数值解,具体形式为:

    \begin{cases} u_t + uu_x - (0.01/\pi)u_{xx} = 0, x\in[-1,1], t\in[0,1]\\u(0,x) = -sin(\pi x), \\ u(t,-1) = u(t,1) = 0\end{cases}

    训练流程

    1.定义网络

            这里定义网络的方法还是和第二篇是一样的,有需要可以去看看第二篇文章

    神经网络学习(二):解常微分方程_lny161224的博客-CSDN博客前言在完成了函数拟合之后,现在考虑常微分方程:给定一个常微分方程,让你求得这个微分方程的近似解(在给定的初值条件下)。以前都是用数学的知识来解决(这让我想起了大一在立人楼上常微分的夜晚),现在有了神经网络,我们用深度学习的方法来解决它试一试。 本次选择的微分方程是训练流程1.定义网络 本来是想使用之前定义的网络,但是觉得以后反正要用就规范化一下这个网络,增加一下可操作性(其实是因为当时我人在办公室笔记本上没有代码就直接重写了),在初始化中增加...https://blog.csdn.net/lny161224/article/details/120499386

            因为输入是二维的(由t,x形成的二维矩阵),所以把输入层改为2即可。

    class Net(nn.Module):
        def __init__(self, NL, NN):
            # NL是有多少层隐藏层
            # NN是每层的神经元数量
            super(Net, self).__init__()
            self.input_layer = nn.Linear(2, NN)
            self.hidden_layer = nn.ModuleList([nn.Linear(NN, NN) for i in range(NL)])
            self.output_layer = nn.Linear(NN, 1)
    
        def forward(self, x):
            o = self.act(self.input_layer(x))
            for i, li in enumerate(self.hidden_layer):
                o = self.act(li(o))
            out = self.output_layer(o)
            return out
    
        def act(self, x):
            return torch.tanh(x)

    2.初始化样本集

            与之前不同的是,这里的样本分为了两类,一类是用于训练用的随机二维样本,另一类是满足初值条件和边界条件的特殊二维样本,所以为了方便整理,定义了一个sample函数来生成样本集

    def sample(size):
        x = torch.cat((torch.rand([size, 1]), torch.full([size, 1], -1) + torch.rand([size, 1]) * 2), dim=1)
        # 因为后面有sin(πx),所以这里得先把x存下来,并且保证和带入初值的x_initial(t,x)的x部分是相同的
        x_init = torch.full([size, 1], -1) + torch.rand([size, 1]) * 2
        x_initial = torch.cat((torch.zeros(size, 1), x_init), dim=1)
        x_boundary_left = torch.cat((torch.rand([size, 1]), torch.full([size, 1], -1)), dim=1)
        x_boundary_right = torch.cat((torch.rand([size, 1]), torch.ones([size, 1])), dim=1)
        return x, x_initial, x_init, x_boundary_left, x_boundary_right

            这里处理上有个手法,因为x\in[-1,1],t\in[0,1]所以,生成t的时候直接用torch.rand()就可以,而对于x,我想到的方法就是从-1开始,生成[0,2]的随机结果然后合在一起,那么结果就是在[-1,1]上的(虽然忘了在哪学的,但是感觉就该这么处理嘻嘻嘻

            这里使用了torch.cat()将两个维度拼在一起,形成了想要的样本集。这里要着重讲一下初值条件,因为初值条件是在t = 0时得到的只跟x有关的函数,所以生成初值样本集的时候一定要注意后面用于计算的u(0,x)sin(\pi x)所对应的x应该是一样的,这里就直接保存了随机结果x_init,这样保证了一致性,且把x_init传了回来,这个在后面会提到为什么会传过来

    3.定义损失

            还是跟(二)一样,将近似函数的损失(误差)和初值、边界条件的误差加在一起作为其损失,处理如下

            loss1 = loss_fn(dt + y_train * dx, (0.01 / np.pi) * dxx)
            loss2 = loss_fn(net(x_initial), -1 * torch.sin(np.pi * x_init))
            loss3 = loss_fn(net(x_left), torch.zeros([size, 1]))
            loss4 = loss_fn(net(x_right), torch.zeros([size, 1]))
            loss = loss1 + loss2 + loss3 + loss4

            这里有两个地方需要注意:

            a)从思路上来说我们应该吧等式两边移到一边,然后使他和0对比作为损失,只不过这样操作的话就需要想loss3loss4一样生成同样size的零矩阵,为了节省空间(实际是为了偷懒),等式两边都放一部分,这样就不需要再额外生成了

            b)注意这个loss2,等式左边是在网络下得到的初值结果,右边是同样的x的情况下得到的值,首先得保证这两个对应的x是一样的,其次需要注意的是右边只需要传入x这一个维度,不需要u(卡了几天没解决都是因为这个地方),因此才需要生成样本的时候传回x_init

    4.函数求偏导

            这里充分利用了autograd自动求梯度来实现,二阶导就是在一阶导的情况下再来一次

            x = Variable(x_train, requires_grad=True)
            d = torch.autograd.grad(net(x), x, grad_outputs=torch.ones_like(net(x)), create_graph=True)
            dt = d[0][:, 0].unsqueeze(-1)
            dx = d[0][:, 1].unsqueeze(-1)
            dxx = torch.autograd.grad(dx, x, grad_outputs=torch.ones_like(dx), create_graph=True)[0][:, 1].unsqueeze(-1)
            # 先把该求的求完然后清零,减少影响(虽然不知道有没有用)
            optimizer.zero_grad()

    5.训练模型

        for i in range(10 ** 4):
            x = Variable(x_train, requires_grad=True)
            d = torch.autograd.grad(net(x), x, grad_outputs=torch.ones_like(net(x)), create_graph=True)
            dt = d[0][:, 0].unsqueeze(-1)
            dx = d[0][:, 1].unsqueeze(-1)
            dxx = torch.autograd.grad(dx, x, grad_outputs=torch.ones_like(dx), create_graph=True)[0][:, 1].unsqueeze(-1)
            # 先把该求的求完然后清零,减少影响(虽然不知道有没有用)
            optimizer.zero_grad()
            y_train = net(x)
            # 为了少生成几个零矩阵,所以计算损失的时候把条件移到两边了
            loss1 = loss_fn(dt + y_train * dx, (0.01 / np.pi) * dxx)
            loss2 = loss_fn(net(x_initial), -1 * torch.sin(np.pi * x_init))
            loss3 = loss_fn(net(x_left), torch.zeros([size, 1]))
            loss4 = loss_fn(net(x_right), torch.zeros([size, 1]))
            loss = loss1 + loss2 + loss3 + loss4
            loss.backward()
            optimizer.step()
            if i % 1000 == 0:
                print(f'step: {i}  loss = {loss.item()}')

    6.数据可视化

            因为自变量成为了2维,所以想要表示u(t,x)得用3D来画图,并把图像保存下来

        T, X = np.meshgrid(t, x, indexing='ij')
        pred_surface = np.reshape(pred, (t.shape[0], x.shape[0]))
        Exact_surface = np.reshape(Exact, (t.shape[0], x.shape[0]))
    
        # plot the approximated values
        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')
        ax.set_zlim([-1, 1])
        ax.plot_surface(T, X, pred_surface, cmap=cm.RdYlBu_r, edgecolor='blue', linewidth=0.0003, antialiased=True)
        ax.set_xlabel('t')
        ax.set_ylabel('x')
        ax.set_zlabel('u')
        plt.savefig('Preddata.png')
        plt.close(fig)

    7.训练结果

            为了知道模型是不是正确的,我找了PINNs里给出的正确解,并且在同样的范围内用模型求解,附图如下:

     (模型与真实结果对比图:左边是真实数据集形成的u(t,x),右边是预测模型生成的u(t,x))

             可以看到效果还是挺好的,些许的不同可能是因为我的样本点选的比较少,可以加大样本点数量达到更好效果

    总结与期望

             利用深度学习来求偏微分方程数值解,只要知道如何将导数表示出来,就可以很好地表示出损失剩下的就只需要交给你的网络,让网络慢慢去学习更新优化参数,最后获得比较满意的结果。

            本文实现了PINNs里Burgers公式求解的代码复现,其实这也是我最开始遇到的问题,前面两篇就是为了解决这个问题做的铺垫(因为模型效果一直不好,所以想到简化问题看能不能实现),结果发现就是因为损失那个地方被sin(πx)整了一手,这里的x并非输入层得到的x(输入层得到的实际是u(t,x),而sin(πx)是输入层第二维的x),总的来说,收获颇丰,很多东西都是自己一点一点摸索出来的,当然也是学习了知乎大佬的解法,慢慢找到正确的解法,再次感谢!

    深度学习求解偏微分方程系列一:Deep Galerkin Method - 知乎我们接下来将用一个系列的文章,介绍使用神经网络的办法求解偏微分方程。这是本系列的第一篇,介绍Deep Galerkin Method (DGM)。我们首先介绍DGM的理论。最后我们使用Python求解一个热传播方程。1. 简介偏微分方…https://zhuanlan.zhihu.com/p/359328643        接下来准备考虑实现的是PINNs中的第二个问题:在数据驱动下的探索发现,即得到了正确的数据集u(t,x),但是偏微分方程中是带了未知参数的情况,通过神经网络来求得参数的方法

    源代码

    import torch
    import torch.nn as nn
    import numpy as np
    import matplotlib.pyplot as plt
    import torch.optim as optim
    from torch.autograd import Variable
    from matplotlib import cm
    import scipy.io
    
    
    class Net(nn.Module):
        def __init__(self, NL, NN):
            # NL是有多少层隐藏层
            # NN是每层的神经元数量
            super(Net, self).__init__()
            self.input_layer = nn.Linear(2, NN)
            self.hidden_layer = nn.ModuleList([nn.Linear(NN, NN) for i in range(NL)])
            self.output_layer = nn.Linear(NN, 1)
    
        def forward(self, x):
            o = self.act(self.input_layer(x))
            for i, li in enumerate(self.hidden_layer):
                o = self.act(li(o))
            out = self.output_layer(o)
            return out
    
        def act(self, x):
            return torch.tanh(x)
    
    
    """
    用来生成拟合函数、初值条件、边界条件的输入
    这里是用torch.rand随机生成,并且根据取值范围调整了表达式
    """
    
    
    def sample(size):
        x = torch.cat((torch.rand([size, 1]), torch.full([size, 1], -1) + torch.rand([size, 1]) * 2), dim=1)
        # 因为后面有sin(πx),所以这里得先把x存下来,并且保证和带入初值的x_initial(t,x)的x部分是相同的
        x_init = torch.full([size, 1], -1) + torch.rand([size, 1]) * 2
        x_initial = torch.cat((torch.zeros(size, 1), x_init), dim=1)
        x_boundary_left = torch.cat((torch.rand([size, 1]), torch.full([size, 1], -1)), dim=1)
        x_boundary_right = torch.cat((torch.rand([size, 1]), torch.ones([size, 1])), dim=1)
        return x, x_initial, x_init, x_boundary_left, x_boundary_right
    
    
    if __name__ == '__main__':
        size = 2000
        lr = 1e-4
        x_train, x_initial, x_init, x_left, x_right = sample(size)
        net = Net(NL=4, NN=20)
        optimizer = optim.Adam(net.parameters(), lr)
        loss_fn = nn.MSELoss(reduction='mean')
        for i in range(10 ** 4):
            x = Variable(x_train, requires_grad=True)
            d = torch.autograd.grad(net(x), x, grad_outputs=torch.ones_like(net(x)), create_graph=True)
            dt = d[0][:, 0].unsqueeze(-1)
            dx = d[0][:, 1].unsqueeze(-1)
            dxx = torch.autograd.grad(dx, x, grad_outputs=torch.ones_like(dx), create_graph=True)[0][:, 1].unsqueeze(-1)
            # 先把该求的求完然后清零,减少影响(虽然不知道有没有用)
            optimizer.zero_grad()
            y_train = net(x)
            # 为了少生成几个零矩阵,所以计算损失的时候把条件移到两边了
            loss1 = loss_fn(dt + y_train * dx, (0.01 / np.pi) * dxx)
            loss2 = loss_fn(net(x_initial), -1 * torch.sin(np.pi * x_init))
            loss3 = loss_fn(net(x_left), torch.zeros([size, 1]))
            loss4 = loss_fn(net(x_right), torch.zeros([size, 1]))
            loss = loss1 + loss2 + loss3 + loss4
            loss.backward()
            optimizer.step()
            if i % 1000 == 0:
                print(f'step: {i}  loss = {loss.item()}')
        """
        接下来引入正确解 画图进行比较
        """
        data = scipy.io.loadmat('burgers_shock.mat')
        t = data['t'].flatten()[:, None]
        x = data['x'].flatten()[:, None]
        Exact = np.real(data['usol']).T
    
        temp = np.empty((2, 1))
        i = 0
        j = 0
        pred = np.zeros((100, 256))
        for _t in t:
            temp[0] = _t
            for _x in x:
                temp[1] = _x
                ctemp = torch.Tensor(temp.reshape(1, -1))
                pred[i][j] = net(ctemp).detach().cpu().numpy()
                j = j + 1
                if j == 256:
                    j = 0
            i = i + 1
        T, X = np.meshgrid(t, x, indexing='ij')
        pred_surface = np.reshape(pred, (t.shape[0], x.shape[0]))
        Exact_surface = np.reshape(Exact, (t.shape[0], x.shape[0]))
    
        # plot the approximated values
        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')
        ax.set_zlim([-1, 1])
        ax.plot_surface(T, X, pred_surface, cmap=cm.RdYlBu_r, edgecolor='blue', linewidth=0.0003, antialiased=True)
        ax.set_xlabel('t')
        ax.set_ylabel('x')
        ax.set_zlabel('u')
        plt.savefig('Preddata.png')
        plt.close(fig)
        # plot the exact solution
        fig = plt.figure()
        ax = fig.add_subplot(projection='3d')
        ax.set_zlim([-1, 1])
        ax.plot_surface(T, X, Exact_surface, cmap=cm.RdYlBu_r, edgecolor='blue', linewidth=0.0003, antialiased=True)
        ax.set_xlabel('t')
        ax.set_ylabel('x')
        ax.set_zlabel('u')
        plt.savefig('Turedata.png')
        plt.close(fig)
    

    展开全文
  • 学习近似算子的深度学习网络可用于一次求解所有相似的偏微分方程,并针对一系列初始条件和边界条件以及物理参数对相同现象进行建模。 1995年由工作表明,浅层网络可以近似操作符算子。由于引入了神经网络,所以此类...

    点上方计算机视觉联盟获取更多干货

    仅作学术分享,不代表本公众号立场,侵权联系删除

    转载于:新智元

    AI博士笔记系列推荐

    周志华《机器学习》手推笔记正式开源!可打印版本附pdf下载链接

    对于特别复杂的偏微分方程,可能需要数百万个CPU小时才能求解出来一个结果。随着问题越来越复杂,从设计更优秀的火箭发动机到模拟气候变化,科学家们需要一个更「聪明」的求解方法。

    随着任务数量的增加,使用当前计算方法来构建通用的日常机器人的成本变得过高,人们正在快速寻求一种解决办法。我们都希望通用机器人可以执行一系列复杂的任务,例如清洁,维护和交付

    你的「高等数学」还好吗?

     

    微分方程是数学中重要的一课。所谓微分方程,就是含有未知函数的导数。一般凡是表示未知函数、未知函数的导数与自变量之间关系的方程,就叫做微分方程。

     

    如果未知函数是一元函数的,就叫做常微分方程;

     

    如果未知函数是多元的,就叫做偏微分方程。

     

    偏微分方程拥有广泛的应用场景,模拟客机在空中的飞行姿势,地震波在地球上的仿真,传染病在人群中扩散的过程,研究基本力和粒子之间的相互作用等场景,工程师、科学家和数学家们都诉诸于偏微分方程来描述涉及许多独立变量的复杂现象。

     

    然而,偏微分方程的求解过程却是异常艰难的,尤其对于计算机来说,只能以最笨拙的方法去求解。

     

    对于特别复杂的偏微分方程,可能需要数百万个CPU小时才能求解出来一个结果。随着问题越来越复杂,从设计更优秀的火箭发动机到模拟气候变化,科学家们需要一个更「聪明」的求解方法。

     

    或许,可以试试神经网络?

     

    最近,研究人员通过神经网络实验证明了它可以比传统的偏微分方程求解器更快地求出近似解。

     

    更牛的是,经过训练的网络,无需再次训练就可以解决一类偏微分方程。

     

    通常,神经网络将数据从一个有限维空间(例如,图像的像素值)映射或转换为另一个有限维空间(例如,将图像分类的数字,例如1代表猫,2代表狗)。

     

    求解偏微分方程的神经网络从无穷大的空间映射到无穷大的空间。

     

     

    偏微分方程的用处和他们的复杂性相伴而生,例如,我们想要观察空气在飞机机翼附近的流动二维透视图,建模人员想知道流体在空间中任何一点(也称为流场)以及在不同时间的速度和压力的话,就需要用到偏微分方程。

     

    考虑到能量、质量和动量守恒定律,特定的偏微分方程,即Navier-Stokes方程可以对这种流体流动进行建模。

     

    在这种情况下,解决方案可能是一个特定公式,可以让开发人员在不同时间计算流场的状态。

     

    偏微分方程常常是很复杂的,以至于无法提供通用的分析解决方案。对于Navier-Stokes方程的最通用形式尤其如此:数学家尚未证明是否存在唯一解,更不用说实际地通过分析找到它们了。

     

    在这些情况下,建模者会转向数值方法,将偏微分方程转换为一组易于处理的代数方程,假定这些方程可保持很小的空间和时间增量。

     

    在超级计算机上,用数值方式解决复杂的偏微分方程可能要花费数月的时间。

     

    而且,如果初始条件或边界条件或所研究系统的几何形状(例如机翼设计)发生了变化,就必须重新开始求解。同样,使用的增量越小(如研究人员所说,网格越细),模型的精度越高,数值求解所需的时间就越长。

     

    神经网络更擅长拟合这样一个黑盒的未知函数,输入是一个向量,而输出是另一个向量。如果存在将一组输入向量映射到一组输出向量的函数,则可以训练网络以学习该映射,两个有限维空间之间的任何函数都可以通过神经网络近似。

     

    2016年,研究人员研究了如何将通常用于图像识别的深度神经网络用于解决偏微分方程。首先,研究人员生成了用于训练网络的数据:一个数值求解器计算了流过xy且大小和方向不同的基本形状(三角形,四边形等)的简单对象上流动的流体的速度场。2D图像编码有关对象几何形状和流体初始条件的信息作为输入,而相应速度场的2D快照作为输出。

     

     

    从无限空间映射到无限空间

     

    相比2016年的工作,这次的研究更有飞跃性的意义,该网络不仅可以学习如何近似函数,还可以学习将函数映射到函数的「运算符」,而且没有「维度爆炸」的困扰。例如,如果其他神经网络或机器学习算法,希望错误率从10%下降到1%,则所需的训练数据量或网络规模可能会成倍爆炸,从而使任务无法实现。

     

    在数学上,操作符的输入输出是没有限制的,例如正弦函数sin(x),输入和输出端是无穷维的,因为x可以是任何值,函数可以是作用于x的任何变换。

     

    学习近似算子的深度学习网络可用于一次求解所有相似的偏微分方程,并针对一系列初始条件和边界条件以及物理参数对相同现象进行建模。

     

     

    1995年由工作表明,浅层网络可以近似操作符算子。由于引入了神经网络,所以此类算子称为神经算子,即实际算子的近似值。

     

    2019年,研究人员提出DeepONet,基于1995年的工作。它的的独特之处在于它的分叉架构,该架构在两个并行网络(一个分支和一个主干)中处理数据。前者在输入端学习一些函数的近似值,而后者在输出端学习相同的函数。

     

    DeepONet将两个网络的输出合并,以了解偏微分方程所需的运算符。训练DeepONet并在每次迭代中调整分支网络和主干网络中的权重,直到整个网络几乎没有出现误差允许范围外的错误为止。

     

     

    DeepONet一旦训练后,就可以模拟操作符,可以在输入端获取代表偏微分方程的数据,其输出为网络训练得到的近似解。

     

    假设您提供了100个样本,这些样本代表了训练数据中没有的初始/边界条件和物理参数,以及需要流场的位置,DeepONet可以在几分之一秒内为您提供流场的状态。

     

    但,DeepONet的训练过程仍然需要消耗大量算力,并且如何提升精确度,以及缩小步长产生更大的计算,也是一个问题。还能更快吗?

    改变观点

    去年,加州理工学院和普渡大学的Anandkumar及其同事建立了一个称为傅立叶神经算子(FNO)的深度神经网络,他们声称这种神经网络具有更快的速度。

     

    他们的网络还将函数映射到函数,从无穷维空间到无穷维空间,并且他们在偏微分方程上测试了它们的神经网络。

     

    他们解决方案的核心是一个傅立叶层。

     

    在他们将训练数据推过神经网络的单层之前,他们先对其进行了傅里叶变换。然后,当图层通过线性运算处理了该数据时,他们将执行傅立叶逆变换,将其转换回原始格式,此转换是著名的傅里叶变换,它将连续函数分解为多个正弦函数。

     

    整个神经网络由几个傅立叶层组成。

     

    事实证明,此过程比DeepONet的计算更直接,并且类似于通过执行称为PDE与某些其他函数之间的卷积的繁琐数学运算来求解PDE。

     

    在傅立叶域中,卷积涉及一个简单的乘法,相当于将经过傅立叶变换的数据通过一层人工神经元(在训练过程中获得的精确权重),然后进行傅立叶逆变换。

     

    因此,最终结果还是FNO学习了整个偏微分方程的操作符,将函数映射到函数。

     

     

    这种方法显著提升了求解速度。

     

    在一个相对简单的示例中,仅需要进行30,000次仿真,就能求出之前提到的Navier-Stokes方程的解,对于每个仿真,FNO花费了几分之一秒的时间,而DeepONet为2.5秒。同样的精度,传统的求解器将花费18个小时。

    数学意义

    两种团队的方法都被证明是成功的,但是与广泛使用神经网络一样,目前尚不清楚它们为什么如此出色以及是否在所有情况下都能如此。Mishra和他的同事现在正在对这两种方法进行全面的数学理解。

     

    经过一年的努力,在2月份,Mishra的团队在Karniadakis的帮助下,对DeepONet架构进行了长达112页的数学分析。他们证明了这种方法是真正通用的,因为它可以将输入端的任何函数集映射到输出端的任何函数集,而不仅仅是PDE,而不必为深入了解Karniadakis定理而做出某些假设网及其1995年的前身。

     

    该团队尚未完成分析FNO的论文,但是Mishra认为它可以比DeepONet更有效地解决某些特定问题。他的团队正在对FNO进行详细的分析,其中包括与DeepONet的比较。

     

    但是,很明显的是,这两种方法都会超越传统的求解器。对于一些无法写出偏微分方程的场景中,神经算子可能是建模此类系统的唯一方法。

     

    这是科学机器学习的未来。

    参考资料:https://www.quantamagazine.org/new-neural-networks-solve-hardest-equations-faster-than-ever-20210419/

    -------------------

    END

    --------------------

    我是王博Kings,985AI博士,华为云专家、CSDN博客专家(人工智能领域优质作者)。单个AI开源项目现在已经获得了2100+标星。现在在做AI相关内容,欢迎一起交流学习、生活各方面的问题,一起加油进步!

    我们微信交流群涵盖以下方向(但并不局限于以下内容):人工智能,计算机视觉,自然语言处理,目标检测,语义分割,自动驾驶,GAN,强化学习,SLAM,人脸检测,最新算法,最新论文,OpenCV,TensorFlow,PyTorch,开源框架,学习方法...

    这是我的私人微信,位置有限,一起进步!

    王博的公众号,欢迎关注,干货多多

    王博Kings的系列手推笔记(附高清PDF下载):

    博士笔记 | 周志华《机器学习》手推笔记第一章思维导图

    博士笔记 | 周志华《机器学习》手推笔记第二章“模型评估与选择”

    博士笔记 | 周志华《机器学习》手推笔记第三章“线性模型”

    博士笔记 | 周志华《机器学习》手推笔记第四章“决策树”

    博士笔记 | 周志华《机器学习》手推笔记第五章“神经网络”

    博士笔记 | 周志华《机器学习》手推笔记第六章支持向量机(上)

    博士笔记 | 周志华《机器学习》手推笔记第六章支持向量机(下)

    博士笔记 | 周志华《机器学习》手推笔记第七章贝叶斯分类(上)

    博士笔记 | 周志华《机器学习》手推笔记第七章贝叶斯分类(下)

    博士笔记 | 周志华《机器学习》手推笔记第八章集成学习(上)

    博士笔记 | 周志华《机器学习》手推笔记第八章集成学习(下)

    博士笔记 | 周志华《机器学习》手推笔记第九章聚类

    博士笔记 | 周志华《机器学习》手推笔记第十章降维与度量学习

    博士笔记 | 周志华《机器学习》手推笔记第十一章稀疏学习

    博士笔记 | 周志华《机器学习》手推笔记第十二章计算学习理论

    博士笔记 | 周志华《机器学习》手推笔记第十三章半监督学习

    博士笔记 | 周志华《机器学习》手推笔记第十四章概率图模型

    点分享

    点收藏

    点点赞

    点在看

    展开全文
  • 转自:新智元来源:quanta,编辑:LRS对于特别复杂的偏微分方程,可能需要数百万个CPU小时才能求解出来一个结果。随着问题越来越复杂,从设计更优秀的火箭发动机到模拟气候变化,科学家们需...
  • 图神经网络解偏微分方程系列(一)1. 标题和概述Learning continuous-time PDEs from sparse(稀疏) data with graph neural networks 使用图神经网络从稀疏数据中学习连续时间偏微分方程这篇文章是使用图神经网络从...
  • 神经网络求解偏微分方程代码分析

    千次阅读 2021-03-20 14:04:15
    下面把我求解偏微分方程的代码分享出来,主要是分享代码思路。这个代码是在求解常微分方程的基础上进行的修改,现在看来有些语句可以换成更高级的表达。若有更好的表达方式欢迎评论/私信! 运行环境:python3.6 + ...
  • 在北京智源大会的“人工智能的数理基础专题论坛”上,北京大学副教授、智源学者董彬做了题为《Learning and Learning to Solve PDEs》的主题演讲(编者注:PDE,Partial Differential Equation,偏微分方程)。...
  • 2.4 一阶隐式微分方程参数表示

    千次阅读 2020-01-16 15:54:55
    文章目录类型一:y=f(x,dydx)y=f(x,\frac{dy}{dx})y=f(x,dxdy​)类型二:x=f(y,dydx)x=f(y,\frac{dy}{dx})x=f(y,dxdy​)一个题解类型三:不显含y的方程F(x,y′)=0F(x,y')=0F...一阶隐式微分方程通常表示:F(x,y,y′...
  • 这是一个二维正交椭圆网格(网​​格)生成器,它通过求解Winslow偏微分方程(Elliptic PDE)来工作。 它能够使用拉伸功能和正交性调整算法来修改网格。 该算法通过使用倾斜的抛物线切线拟合器计算曲线斜率来工作。 ...
  • 神经网络求解二阶常微分方程(代码)

    千次阅读 2020-11-29 20:24:27
    相关的理论推导请参考: 神经网络求解二阶常微分方程. 以下是程序分享 import os os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' import tensorflow as tf import matplotlib.pyplot as plt import numpy as np import ...
  • 微分方程建模详解

    千次阅读 2021-01-25 00:00:00
    微分方程建模
  • 神经网络求解二阶常微分方程

    千次阅读 热门讨论 2020-11-20 15:10:41
    需要让我使用深度神经网络求解偏微分方程。在相关调研过程中,CSDN上作者Trytobenice分享过相关的程序源码。基于相关程序源码,我将他的一阶常微分方程求解扩充到二阶常微分方程求解。并且按照此方法可以求解高阶常...
  • Python的Scipy库解微分方程

    千次阅读 2020-01-16 01:08:02
    微分方程: 初始值: 问题: 求解其他三个参数: 代码实现: import numpy as np from numpy import zeros, linspace, arange from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_...
  • 点击上方,选择星标或置顶,不定期资源大放送!阅读大概需要15分钟Follow小博主,每天更新前沿干货转载自:量子位AI也能解方程了?是的,它们不仅能解方程,还能“找到”方程!今天我们就简单...
  • 本文以一个微分方程为例求解该微分方程的数值解,并给出模拟数值曲线。普通的微分方程(组)亦可使用此方法求解并作图。在本文后面部分,给出了python内置函数求解,并比较了二者的区别。
  • 在完成了函数拟合之后,现在考虑常微分方程:给定一个常微分方程,让你求得这个微分方程的近似解(在给定的初值条件下)。以前都是用数学的知识来解决(这让我想起了大一在立人楼上常微分的夜晚),现在有了神经网络...
  • 在微分方程中,对一个变量微分称为常微分方程,对多个变量微分称为偏微分方程。 劳伦兹(Lorenz)方程式是由大气方程简化而来的,方程式如下: 其中,a,b,c为特定参数,在某些特定数值中可以观察到一些现象。下面...
  • 通过Julia语言,实现常用的非线性方程和常微分方程的求解方法,包括二分法,Newton法,Runge-Kutta法(龙格-库塔法),Euler法(欧拉法)等。
  • 前面几篇博客介绍了神经网络应用到积分、一元N阶微分的原理、方法并实践了可行性,取得了较好的拟合效果,现在针对偏微分方程PDE进行最后的攻关,完成该部分攻关后即基本掌握了神经网络应用到方程求解的原理方法以及...
  • 书接上回,我们在这里讨论一下微分方程模型,也是预测模型的最后一节,以后有想到的再补上、()拟合优度对于非线性情况已经没有意义了。。 分类 ...多变量:偏微分方程,因此,微分方程指的是又有
  • 本文研究基于线性常微分方程的福建省艾滋病传播数学模型的建立和预测分析问题。对疾病传播中传统的SI模型进行建模和分析,提出SEID模型,新增和文化程度相关的日接触率参数、和当地卫生状况、医疗水平相关的发病率、...
  • 本博客对更复杂的二阶微分问题进行神经网络求解,问题示例参考博客 问题描述 定义域选取[0,2] 模型代码 神经网络模拟φ(x),利用自动微分得到二阶、一阶微分,代入表达式后作为loss进行训练即可。该方法适用于N阶...
  • 第三场微分方程建模讲座笔记 主讲人:董程栋,上海财经大学数学学院 微分方程: 定义:联系着自变量,未知函数与它的导数之间的关系式 物体冷却过程中的数学模型 牛顿冷却定律:物体温度变化速度与物体和介质...
  • 数学建模(三)模型拟合

    千次阅读 2019-08-14 22:06:33
    在数学建模过程中,需要根据不同目的分析数据,当问题很复杂难以建立能够解释该特殊情形的模型时,如果子模型涉及偏微分方程,并且没有封闭解的时候,那么再以此构造一个主模型时将很难得到解,这个时候就需要进行...
  • MATLAB解微分方程
  • 0.导语Scipy是一个用于数学、科学、工程领域的常用软件包,可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵,使Numpy和Scipy协同工作,高效解决问题。Scipy...
  • 微分方程建模 微分方程是数学建模的重要方法,因为许多实际问题的数学描述将导致求解微分方程的定解问题。 把实际问题化为求解微分方程的定解问题,大体上可以分为以下步骤 根据实际要求确定要研究的量(自变量,...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,452
精华内容 980
关键字:

偏微分方程参数拟合