精华内容
下载资源
问答
  • 隐式求解偏微分方程
    千次阅读
    2021-01-30 02:15:19

    我有一个偏微分方程组(PDE),特别是应用于传热和对流的扩散-平流-反应方程,我用有限差分法求解。

    在我的模拟环境中,我有很多不同的部分,比如管道,储能器,热交换器等等。。。根据零件的不同,每个零件的PDE可以是1D(在大多数情况下)或2D(约占零件的5%)或3D(很少)。由于这些不同的形状,为整个系统构造一个三对角(或五角、间隔-…)太复杂了。

    这个问题发表在SO上而不是计算科学上,因为它是关于如何使用pde的求解算法的,而不是关于pde的求解方法。在

    因此,目前每个部件都有一个diff函数,它返回给定输入温度的微分。由于材料特性等与温度(和流量)有关,因此PDE是非线性的,但由于系数滞后于系数而被认为是线性的(在每个步骤使用起始温度计算温度和流量相关值,对于在每次迭代中重新计算它们的隐式方法)。diff函数返回的微分的形状与零件的值数组形状相同,例如,对于具有10个网格点的一维管道,这将导致形状(10, )的差分数组。扩散系数和其他系数被认为是diff函数的内部参数,对解算器不可用。因此解算器将只知道差分、当前零件温度和步长。

    有没有什么方法可以用python一次只执行一个步骤,使用一个专门用于求解pde的算法来求解这些pde?(以及优选为scipy/numpy的一部分并且更优选地已经由numba支持的算法。)

    到目前为止我考虑过的事情:scipy.integrate:仅对ode可行,而PDE可能不在其中。对吗?至少我不能用它来做PDE的好结果。除此之外,scipy.integrate并不是只用于计算一个步骤的。在

    np.linalg.solve:似乎是最好的方法,如果你有三对角矩阵和线性方程组。但由于我既没有三对角矩阵,也没有线性方程(线性化了,但系数滞后意味着必须在迭代过程中更新它们)。有没有办法在不增加计算成本的情况下仍然使用这个方法?在

    scipy.optimize:可行且易于使用。下面我的最小工作示例将使用此集成一个解决方案。但是这个是否总是适用于pde并且这通常是隐式解决偏微分方程的正确方法,还是仅仅在这个例子中起作用而我走错了方向?

    一个简单迭代函数的自身积分。也显示在下面的例子中。与scipy.optimize相同:在这个简单的例子中,这是有效的,还是对pde一般都有效?在

    最小工作示例import time

    from scipy import optimize

    def diff(T): # simply example differential function

    time.sleep(0.0001) # dummy timer to slow down the calculation

    return np.sqrt(T)

    def euler_bdf(y, yprev, h):

    return (y - yprev - h*diff(y))

    def crank_nicolson(y, yprev, h): # diff(yprev) will be outsourced for performance

    return (y - yprev - h / 2 * (diff(y) + diff(yprev)))

    def iterate_eulerbdf(y0, h, rtol, diff):

    y = y0

    f = np.ones_like(y)

    i = 0

    while np.any(f > rtol):

    y_old = y

    y = y0 + h * diff(y)

    f = (y / y_old - 1) # this will be replaced by least mean squares

    i += 1

    return y, i

    def iterate_cranknicolson(y0, h, rtol, diff):

    y = y0

    diff_old = diff(y)

    f = np.ones_like(y)

    i = 0

    while np.any(f > rtol):

    y_old = y

    y = y0 + h / 2 * (diff_old + diff(y))

    f = (y / y_old - 1)

    i += 1

    return y, i

    T0 = np.arange(10, 20).astype(float) # starting temperature

    eulerbdf_hybr = optimize.root(euler_bdf, T0, args=(T0, 1.), method='hybr', tol=1e-6)

    eulerbdf_dfsane = optimize.root(euler_bdf, T0, args=(T0, 1.), method='df-sane', tol=1e-6)

    eulerbdf_it = iterate_eulerbdf(T0, 1., 1e-6, diff)

    cn_hybr = optimize.root(crank_nicolson, T0, args=(T0, 1.), method='hybr', tol=1e-6)

    cn_dfsane = optimize.root(crank_nicolson, T0, args=(T0, 1.), method='df-sane', tol=1e-6)

    cn_it = iterate_cranknicolson(T0, 1., 1e-6, diff)

    TL,DR问题摘要

    所有的结果似乎都是好的,但这种方法对一般的实偏微分方程行吗?

    推荐的方法是什么?

    对于各种不同形状和尺寸的pde,从(5,)到超过(100, 100, 100),哪种方法最有效?因为scipy.integrate目前{}似乎是最快的,但我怀疑这是否适用于更大的问题。

    尤其是:还有更好的方法吗?

    提前谢谢!在

    更多相关内容
  • 2、古典隐式格式求解抛物型偏微分方程(一维热传导方程) 3、Crank-Nicolson隐式格式求解抛物型偏微分方程 4、正方形区域Laplace方程Diriclet问题的求解 如: function [U x t]=PDEParabolicClassicalExplicit(uX,uT...
  • 使用差分法求解一维的热传导偏微分方程,也可以求解类似热传导的偏微分方程。分别推导了显式和隐式的差分离散格式,并使用matlab编写了显式和隐式求解代码。压缩包中包含了详细的推导文档和带有注释的代码,可供...
  • 【内容介绍】本资源主要利用MATLAB的实时脚本编程实现了抛物型偏微分方程数值求解,以图-文-代码三者互相嵌套的形式介绍实现过程,一目了然。包括对迭代的误差分析。 【适用对象】工科生、数学专业等。 【算法涵盖】...
  • 为研究求解微分方程的近似解问题,采用理论分析和实例分析的方法,将常微分方程的求近似解问题转化为遗传算法的函数优化问题,借助Matlab遗传算法工具箱实现对常微分方程求解,并以室内温度摆动问题进行实例分析....
  • 文章目录(1)偏微分方程的类型(二阶)(2)抛物线型1.显式法2.Crank-Nicholson隐式算法 (3)双曲线型(4)椭圆型 (1)偏微分方程的类型(二阶) a∂2u∂x2+b∂2u∂y∂x+c∂2u∂x2+d∂u∂x+e∂u∂y+fu+g=0a\frac{\partial^2u}{\...

    (1)偏微分方程的类型(二阶)

    a ∂ 2 u ∂ x 2 + b ∂ 2 u ∂ y ∂ x + c ∂ 2 u ∂ x 2 + d ∂ u ∂ x + e ∂ u ∂ y + f u + g = 0 a\frac{\partial^2u}{\partial x^2}+b\frac{\partial^2u}{\partial y\partial x}+c\frac{\partial^2u}{\partial x^2}+d\frac{\partial u}{\partial x}+e\frac{\partial u}{\partial y}+fu+g=0 ax22u+byx2u+cx22u+dxu+eyu+fu+g=0

    • b 2 − 4 a c < 0 b^2-4ac<0 b24ac<0 椭圆
    • b 2 − 4 a c = 0 b^2-4ac=0 b24ac=0 抛物线
    • b 2 − 4 a c > 0 b^2-4ac>0 b24ac>0 双曲线

    (2)抛物线型

    1.显式法

    • 求解思想:通过差分的方法一排一排向上推。
    • 做划分并代入方程 u i , j + 1 − u i , j k = u i + 1 , j − 2 u i , j + u i − 1 , j h 2    ( Δ x = h , Δ t = k ) \frac{u_{i,j+1}-u_{i,j}}{k}=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}~~(\Delta x=h,\Delta t=k) kui,j+1ui,j=h2ui+1,j2ui,j+ui1,j  (Δx=h,Δt=k)
    • 通过化简得到 u i , j + 1 = r u i − 1 , j + ( 1 − 2 r ) u i , j + r u i + 1 , j    ( r = k h 2 ) u_{i,j+1}=ru_{i-1,j}+(1-2r)u_{i,j}+ru_{i+1,j}~~(r=\frac{k}{h^2}) ui,j+1=rui1,j+(12r)ui,j+rui+1,j  (r=h2k)
    • 具体推的步骤大概如下:
      • 由于已知 u ( x , 0 ) = f ( x ) u(x,0)=f(x) u(x,0)=f(x),因此相当于知道 u 0 , 0 , u 1 , 0 , u 2 , 0 … u_{0,0},u_{1,0},u_{2,0}\dots u0,0,u1,0,u2,0
      • 通过上面的公式就可以推出来 u 1 , 1 , u 2 , 1 , u 3 , 1 … u_{1,1},u_{2,1},u_{3,1}\dots u1,1,u2,1,u3,1,注意由于已知左边界和右边界,因此 u 0 , 1 u_{0,1} u0,1也知道,所以第二排就可以全部推出来。
      • 通过上面的方式可以求出区域内全部的数值解。

    2.Crank-Nicholson隐式算法

    • 求解思想:也是一排一排向上推,但是这次是使用线性方程组一次性求出一排。
    • 这里采用相同的划分方式,但是代入不同的差分方程
    • 通过化简得到
    • 具体推的步骤大概如下:
      • 由于已知 u ( x , 0 ) = f ( x ) u(x,0)=f(x) u(x,0)=f(x),因此相当于知道 u 0 , 0 , u 1 , 0 , u 2 , 0 … u_{0,0},u_{1,0},u_{2,0}\dots u0,0,u1,0,u2,0
      • 通过上面的公式就可以推出来方程 − r u i − 1 , 1 + ( 2 + 2 r ) u i , 1 − r u i + 1 , 1 = c     ( c 是 一 个 常 数 ) -ru_{i-1,1}+(2+2r)u_{i,1}-ru_{i+1,1}=c~~~(c是一个常数) rui1,1+(2+2r)ui,1rui+1,1=c   (c),注意由于已知左边界和右边界,所以这个其实就转化成在一维上的差分问题,最后列出全部的方程构成方程组求解即可。
      • 通过上面的方式可以求出区域内全部的数值解。
    • 一个例子:
    • 做划分并且代入差分方程
      k = 0.01 , h = 0.1 k=0.01,h=0.1 k=0.01,h=0.1
      − u i − 1 , j + 1 + 4 u i , j + 1 − u i + 1 , j + 1 = u i − 1 , j + u i + 1 , j -u_{i-1,j+1}+4u_{i,j+1}-u_{i+1,j+1}=u_{i-1,j}+u_{i+1,j} ui1,j+1+4ui,j+1ui+1,j+1=ui1,j+ui+1,j
    • 进行求解(这里利用了对称性,在 x = 0.5 x=0.5 x=0.5 两边是对称的,将 j = 0 j=0 j=0隐去,并根据对称性将 u 6 u_6 u6替换成 u 4 u_4 u4)

    (3)双曲线型

    • 得到的差分方程为 ( r = k h r=\frac{k}{h} r=hk注意和之前的定义不同):
    • 划分需要满足一定的条件 k h ≤ 1 c \frac{k}{h}\le\frac{1}c{} hkc1
    • 具体求解按照之前类似的方法即可。

    (4)椭圆型

    • 得到的差分方程为 (这里取 k k k h h h相等):
    • 求解
      • 可以采用类似之前的隐式或者显式方法求解。
      • 可以采用迭代法求解,比如雅克比迭代,转换成下面的迭代式
    展开全文
  • Python小白的数学建模课-11.偏微分方程数值解法

    千次阅读 多人点赞 2021-08-17 14:10:06
    本文采用有限差分法求解偏微分方程,通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程等常见类型的偏微分方程的数值解法,给出了全部例程和运行结果。。 欢迎关注『Python小白...

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

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

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

    欢迎关注『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数模笔记-模拟退火算法

    展开全文
  • 利用MATLAB求解偏微分方程的有限差分法算例,可自行编写方程求解,只需更换方程即可。
  • 差分算法(求解偏微分方程) 定义 差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分,将自变量的连续变化 区域用有限离散点...

    差分算法(求解偏微分方程)

    差分算法是数学建模比赛中的一种十分常见的代码,在2018A题和2020A中均用到一维热传导模型,模型的求解用的就是差分算法,具体如何解可以自己去查看相关论文。

    定义

    差分方法又称为有限差分方法或网格法,是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是:先对求解区域作网格剖分,将自变量的连续变化 区域用有限离散点(网格点)集代替;将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替;通过用网格点上函数的差商代替导数,将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组(称为差分格式)。如果差分格式有 解,且当网格无限变小时其解收敛于原微分方程定解问题的解,则差分格式的解就作为 原问题的近似解(数值解)。因此,用差分方法求偏微分方程定解问题一般需要解决一下问题:

    1. 选取网络;
    2. 对微分方程及定解条件选择差分近似,列出差分格式;
    3. 求解差分格式;
    4. 讨论差分格式解对于微分方程解的收敛性及误差估计

    算法详解

    请添加图片描述

    因此,只要确定了步长,我们就可以将连续变化的自变量用有限离散点来表示

    请添加图片描述

    对于(9-3)的式子,为了方便计算,我们用差分来表示偏微分方程
    ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = f ( x , y ) \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y) x22u+y22u=f(x,y)
    由于式子中进行了两次偏导,因此我们一步一步进行分析。

    进行一次差分,我们用向前差分
    ∂ u ∂ x = u ( k + 1 , j ) − u ( k , j ) h \frac{\partial u}{\partial x}=\frac{u(k+1,j)-u(k,j)}{h} xu=hu(k+1,j)u(k,j)
    由于向前差分有误差,如果我们进行两次向前差分的话,计算的误差可能会增大,因此,第二次偏导我们选择向后差分。即我们混合向前差分、向后差分来近似代替两次偏导。

    因此,第二次我们用向后差分
    ∂ 2 u ∂ x 2 = ∂ ∂ x ( u ( k + 1 , j ) − u ( k , j ) h ) = ∂ ∂ x ( u ( k + 1 , j ) h ) − ∂ ∂ x ( u ( k , j ) h ) \begin{aligned} \frac{\partial^2u}{\partial x^2}&=\frac{\partial}{\partial x}(\frac{u(k+1,j)-u(k,j)}{h})\\ & = \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})-\frac{\partial}{\partial x}(\frac{u(k,j)}{h}) \end{aligned} x22u=x(hu(k+1,j)u(k,j))=x(hu(k+1,j))x(hu(k,j))

    ∂ ∂ x ( u ( k + 1 , j ) h ) = u ( k + 1 , j ) − u ( k , j ) h ⋅ h − 0 h 2 = u ( k + 1 , j ) − u ( k , j ) h 2 \begin{aligned} \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})&=\frac{\frac{u(k+1,j)-u(k,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k+1,j)-u(k,j)}{h^2} \end{aligned} x(hu(k+1,j))=h2hu(k+1,j)u(k,j)h0=h2u(k+1,j)u(k,j)

    ∂ ∂ x ( u ( k , j ) h ) = u ( k , j ) − u ( k − 1 , j ) h ⋅ h − 0 h 2 = u ( k , j ) − u ( k − 1 , j ) h 2 \begin{aligned} \frac{\partial}{\partial x}(\frac{u(k,j)}{h})&=\frac{\frac{u(k,j)-u(k-1,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k,j)-u(k-1,j)}{h^2} \end{aligned} x(hu(k,j))=h2hu(k,j)u(k1,j)h0=h2u(k,j)u(k1,j)
    综上,
    ∂ 2 u ∂ x 2 = ∂ ∂ x ( u ( k + 1 , j ) − u ( k , j ) h ) = ∂ ∂ x ( u ( k + 1 , j ) h ) − ∂ ∂ x ( u ( k , j ) h ) = u ( k + 1 , j ) − u ( k , j ) h 2 − u ( k , j ) − u ( k − 1 , j ) h 2 = u ( k + 1 , j ) − 2 u ( k , j ) + u ( k − 1 , j ) h 2 \begin{aligned} \frac{\partial^2u}{\partial x^2}&=\frac{\partial}{\partial x}(\frac{u(k+1,j)-u(k,j)}{h})\\ & = \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})-\frac{\partial}{\partial x}(\frac{u(k,j)}{h})\\ &=\frac{u(k+1,j)-u(k,j)}{h^2}-\frac{u(k,j)-u(k-1,j)}{h^2}\\ &= \frac{u(k+1,j)-2u(k,j)+u(k-1,j)}{h^2} \end{aligned} x22u=x(hu(k+1,j)u(k,j))=x(hu(k+1,j))x(hu(k,j))=h2u(k+1,j)u(k,j)h2u(k,j)u(k1,j)=h2u(k+1,j)2u(k,j)+u(k1,j)
    同理可得
    ∂ 2 u ∂ y 2 = u ( k , j + 1 ) − 2 u ( k , j ) + u ( k , j − 1 ) τ 2 \frac{\partial^2u}{\partial y^2}=\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2} y22u=τ2u(k,j+1)2u(k,j)+u(k,j1)
    所以原方程就变为
    ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = f ( x , y ) u ( k + 1 , j ) − 2 u ( k , j ) + u ( k − 1 , j ) h 2 + u ( k , j + 1 ) − 2 u ( k , j ) + u ( k , j − 1 ) τ 2 = f ( x , y ) \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y)\\ \frac{u(k+1,j)-2u(k,j)+u(k-1,j)}{h^2}+\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2}=f(x,y) x22u+y22u=f(x,y)h2u(k+1,j)2u(k,j)+u(k1,j)+τ2u(k,j+1)2u(k,j)+u(k,j1)=f(x,y)
    以上的进行通过混合向前差分、向后差分的方法求二阶偏导的方法实际上就是二阶中心差分

    请添加图片描述
    请添加图片描述

    常见的差分

    向前差分

    函数的前向差分通常简称为函数的差分。对于函数f(x),如果在等距节点:
    x k = x 0 + k h ( k = 0 , 1 , ⋯   , n ) x_k=x_0+kh(k=0,1,\cdots,n) xk=x0+kh(k=0,1,,n)

    Δ f ( x k ) = f ( x k + 1 ) − f ( x k ) \Delta f(x_k)=f(x_{k+1})-f(x_k) Δf(xk)=f(xk+1)f(xk)

    则称Δf(x),函数在每个小区间上的增量 y ( k + 1 ) − y ( k ) y(k+1)-y(k) y(k+1)y(k)为f(x)的一阶前向差分。在微积分学中的有限差分(finite differences),前向差分通常是微分在离散的函数中的等效运算。差分方程的解法也与微分方程的解法相似。当是多项式时,前向差分为Delta算子,一种线性算子。前向差分会将多项式阶数降低1
    由于向前差分可以从已知的值求解出结果,所以我们称向前差分为显式差分。

    向后差分

    向后差分为隐式差分方法,无条件稳定。
    对于函数 f ( x k ) f(x_k) f(xk),一阶向后差分为:
    Δ f ( x k ) = f ( x k ) − f ( x k − 1 ) \Delta f(x_k)=f(x_k)-f(x_{k-1}) Δf(xk)=f(xk)f(xk1)

    Crank-Nicolson差分 (CN差分)

    Crank-Nicolson差分格式又称为中心差分格式。Crank-Nicolson方法式显式方法和隐式方法的结合,式无条件稳定的方法,公式看起来复杂,但是考虑到提高的精度和保证的稳定性。
    Δ f ( x k ) = 1 2 ( f ( x k + 1 ) − f ( x k ) + f ( x k ) − f ( x k − 1 ) ) = 1 2 ( f ( x k + 1 ) − f ( x k − 1 ) ) \begin{aligned} \Delta f(x_k)&=\frac{1}{2}(f(x_{k+1})-f(x_k)+f(x_k)-f(x_{k-1}))\\ &=\frac{1}{2}(f(x_{k+1})-f(x_{k-1})) \end{aligned} Δf(xk)=21(f(xk+1)f(xk)+f(xk)f(xk1))=21(f(xk+1)f(xk1))
    对于扩散方程(包括许多其他方程),可以证明Crank-Nicolson方法无条件稳定。但是,如果时间步长与空间步长平方的比值过大(一般地,大于1/2),近似解中将存在虚假的振荡或衰减。基于这个原因,当要求大时间步或高空间分辨率的时候,往往会采用数值精确较差的后向欧拉方法进行计算,这样即可以保证稳定,又避免了解的伪振荡。

    例题

    用五点菱形个格式求解Laplace方程第一边值问题
    { ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 , ( x , y ) ∈ Ω u ( x , y ) ∣ ( x , y ) ∈ Γ = lg ⁡ [ ( 1 + x ) 2 + y 2 ] , Γ = ∂ Ω , \begin{cases} \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=0,&(x,y)\in \Omega\\ u(x,y)|_{(x,y)\in \Gamma}=\lg [(1+x)^2+y^2],&\Gamma=\partial \Omega, \end{cases} {x22u+y22u=0,u(x,y)(x,y)Γ=lg[(1+x)2+y2],(x,y)ΩΓ=Ω,
    其中 Ω = { ( x , y ) ∣ 0 ≤ x , y ≤ 1 } \Omega=\left\{ (x,y)|0 \leq x,y \leq 1 \right\} Ω={(x,y)0x,y1},取 h = τ = 1 3 h=\tau=\frac{1}{3} h=τ=31

    解:

    请添加图片描述

    根据题意,我们画出网格出来,正好构成四个五点菱形,即得到四个方程,我们将线性方程组写出来
    { 1 h 2 ( u 1 , 2 + u 2 , 1 + u 1 , 0 + u 0 , 1 − 4 u ( 1 , 1 ) ) = 0 1 h 2 ( u 2 , 2 + u 3 , 1 + u 2 , 0 + u 1 , 1 − 4 u ( 2 , 1 ) ) = 0 1 h 2 ( u 1 , 3 + u 2 , 2 + u 1 , 1 + u 0 , 2 − 4 u ( 1 , 2 ) ) = 0 1 h 2 ( u 2 , 3 + u 3 , 2 + u 2 , 1 + u 1 , 2 − 4 u ( 2 , 2 ) ) = 0 \begin{cases} \frac{1}{h^2}(u_{1,2}+u_{2,1}+u_{1,0}+u_{0,1}-4u(1,1))=0\\ \frac{1}{h^2}(u_{2,2}+u_{3,1}+u_{2,0}+u_{1,1}-4u(2,1))=0\\ \frac{1}{h^2}(u_{1,3}+u_{2,2}+u_{1,1}+u_{0,2}-4u(1,2))=0\\ \frac{1}{h^2}(u_{2,3}+u_{3,2}+u_{2,1}+u_{1,2}-4u(2,2))=0 \end{cases} h21(u1,2+u2,1+u1,0+u0,14u(1,1))=0h21(u2,2+u3,1+u2,0+u1,14u(2,1))=0h21(u1,3+u2,2+u1,1+u0,24u(1,2))=0h21(u2,3+u3,2+u2,1+u1,24u(2,2))=0
    线性方程组又可以化简成
    [ − 4 1 1 0 1 − 4 0 1 1 0 − 4 1 0 1 1 − 4 ] [ u 1 , 1 u 1 , 2 u 2 , 1 u 2 , 2 ] = − [ u 1 , 0 + u 0 , 1 u 3 , 1 + u 2 , 0 u 1 , 3 + u 0 , 2 u 2 , 3 + u 3 , 2 ] \begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} 4110140110410114u1,1u1,2u2,1u2,2=u1,0+u0,1u3,1+u2,0u1,3+u0,2u2,3+u3,2

    解非齐次线性方程组求得:
    u 1 = 0.6348 , u 2 = 1.06 , u 3 = 0.7985 , u 4 = 1.1698 u_1=0.6348, u_2=1.06, u_3=0.7985, u_4=1.1698 u1=0.6348,u2=1.06,u3=0.7985,u4=1.1698
    计算的MATLAB程序如下

    clc;
    clear;
    f1 = @(x)2 * log(1+x);
    f2 = @(x)log((1+x).^2+1);
    f3 = @(y)log(1+y.^2);
    f4=  @(y)log(4+y.^2);
    u=zeros(4);
    m=4;% 总列数
    n=4;% 总行数
    h=1/3;
    u(1,1:m)=feval(f3,0:h:(m-1)*h)';
    u(n,1:m)=feval(f4,0:h:(m-1)*h)';
    u(1:n,1)=feval(f1,0:h:(n-1)*h);
    u(1:n,m)=feval(f2,0:h:(n-1)*h);
    b = -[u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3)];
    a = [-4,1,1,0;1,-4,0,1;1,0,-4,1;0,1,1,-4];
    x=a\b;
    

    代码解读

    观察线性方程:
    [ − 4 1 1 0 1 − 4 0 1 1 0 − 4 1 0 1 1 − 4 ] [ u 1 , 1 u 1 , 2 u 2 , 1 u 2 , 2 ] = − [ u 1 , 0 + u 0 , 1 u 3 , 1 + u 2 , 0 u 1 , 3 + u 0 , 2 u 2 , 3 + u 3 , 2 ] \begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} 4110140110410114u1,1u1,2u2,1u2,2=u1,0+u0,1u3,1+u2,0u1,3+u0,2u2,3+u3,2
    这个形式为最小二乘法的形式

    参考Matlab中的线性最小二乘的标准型:
    min ⁡ A ∥ R A − Y ∥ 2 2 \min _A \Vert RA-Y \Vert_2^2 AminRAY22
    所以我们只需要求出等式右边的式子,那么我们就可以解出等式中间的向量

    现在我们来解等式右边的向量,观察可知,为定义域的边界

    注意这里从1开始,与上式子有些不同
    [ u 1 , 1 u 1 , 2 u 1 , 3 u 1 , 4 u 2 , 1 u 2 , 2 u 2 , 3 u 2 , 4 u 3 , 1 u 3 , 2 u 3 , 3 u 3 , 4 u 4 , 1 u 4 , 2 u 4 , 3 u 4 , 4 ] \begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\ u_{2,1} & u_{2,2} & u_{2,3} & u_{2,4}\\ u_{3,1} & u_{3,2} & u_{3,3} & u_{3,4}\\ u_{4,1} & u_{4,2} & u_{4,3} & u_{4,4} \end{bmatrix} u1,1u2,1u3,1u4,1u1,2u2,2u3,2u4,2u1,3u2,3u3,3u4,3u1,4u2,4u3,4u4,4
    因此我们只需要解出边界的值,就可以把b表示出来


    [ u 1 , 1 u 1 , 2 u 1 , 3 u 1 , 4 u 2 , 1 u 2 , 2 u 2 , 3 u 2 , 4 u 3 , 1 u 3 , 2 u 3 , 3 u 3 , 4 u 4 , 1 u 4 , 2 u 4 , 3 u 4 , 4 ] \begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\ u_{2,1} & u_{2,2} & u_{2,3} & u_{2,4}\\ u_{3,1} & u_{3,2} & u_{3,3} & u_{3,4}\\ u_{4,1} & u_{4,2} & u_{4,3} & u_{4,4} \end{bmatrix} u1,1u2,1u3,1u4,1u1,2u2,2u3,2u4,2u1,3u2,3u3,3u4,3u1,4u2,4u3,4u4,4
    因此我们只需要解出边界的值,就可以把b表示出来,然后用最小二乘法解出中间的变量,那么每个点值都接出来了。

    参考文献:数值分析(原书第2版)_Timothy Sauer_机械工业出版社

    展开全文
  • 偏微分方程数值求解 -- ING

    千次阅读 2019-07-06 15:53:35
    偏微分方程的数值解法,本质是偏微分方程的近似求解。虽然是近似,但整个近似过程满足严格的数学要求,如稳定性、收敛性、相容性等,并非随意的近似。由于时间,精力以及个人需求,仅秉着“理解不深究”的态度,进行...
  • 5、显式法、Crank-Nicholson隐式法(抛物型偏微分方程) [图片] 上图疑似有误,应为 Crank-Nicholson 隐式法 边值为u 显式法 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1KbPzh15-...
  • 本文参考的是来自mooc上北京师范大学彭芳麟老师的计算物理基础基础知识偏微分方程的三种类型椭圆型 初始条件:无抛物型 初始条件:初始温度分布双曲型初始条件:初始位移与初始速度边界条件Dirchlet边界条件区域边界...
  • 差分法数值求解热传导偏微分方程

    千次阅读 2022-03-03 21:47:47
    使用差分法求解一维的热传导偏微分方程,也可以求解类似热传导的偏微分方程。分别推导了显式和隐式的差分离散格式,并使用matlab编写了显式和隐式求解代码。压缩包中包含了详细的推导文档和带有注释的代码,可供...
  • 先求通解再确定特解,是求常微分方程定解问题采用的方法,都某些偏微分方程,也能通过积分求出通解,进而确定出满足定解条件的特解。 两个自变量的一阶线性偏微分方程 今有两个自变量的一阶线性偏微分方程。 a(x,y)...
  • Matlab求解微分方程(组)及偏微分方程(组) 共12页.pdf
  • 来源:新浪了凡春秋的博客在科学技术各领域中,有很多问题都可以归结为偏微分方程问题。...下面的几个简单例子,将为大家介绍如何利用Matlab中的PDE工具箱进行偏微分方程求解!抛物线型受热金...
  • matlab:使用欧拉方法求解微分方程

    千次阅读 2022-04-18 10:27:59
    %欧拉方法求解微分方程 function [t,y] = my_euler(f, t0, tf , y0, h) %f-函数; t0,tf:区间; y0,初值;h 步长 M = floor((tf - t0)/h); % 离散点的个数 t = zeros(M +1, 1); y = zeros(M +1, 1); t =...
  • 里面的都能用Pie Eais ae h t lat BLep话k气飞“目基回古班比比,一生热件身的些懂1身定回量1古典显式格式稳定情况于吉氏,一生售挥导方栏的解里点丁空真ecin.com2、古典隐式格式求解抛物型偏微分方程(一维热传导方程)...
  • 偏微分方程数值解法MATLAB源码》由会员分享,可在线阅读,更多相关《偏微分方程数值解法MATLAB源码(27页珍藏版)》请在人人文库网上搜索。1、源码【更新完毕】偏微分方程数值解法的MATLAB原创 说明:由于偏微分的...
  • 数学建模入门-matlab实现偏微分方程数值解

    万次阅读 多人点赞 2019-03-17 22:23:36
    文章目录前言调用示例例题求解命令介绍具体实现步骤1:化标准式步骤 2:编写偏微分方程的系数向量函数步骤3:编写初始条件函数步骤 4:编写边界条件函数步骤 5: 取点主程序 前言 在python3安装fipy失败之后,懒得...
  • 包含椭圆,抛物线,双曲线偏微分方程数值解法,隐式格式,显示格式等,应用于大学偏微分方程数值解报告的撰写
  • 一般地,n个自变量的二阶线性偏微分方程可表示为 ∑i,j=1naij(x1,⋅⋅⋅,xn)∂2u∂xi∂xj+∑j=1nbj(x1,⋅⋅⋅,xn)∂u∂xj+c(x1,⋅⋅⋅,xn)u=f(x1,⋅⋅⋅,xn)(1) \sum_{i,j=1}^na_{ij}(x_1,···,x_n)\frac{\partial...
  • MATLAB源码--古典显式格式求解抛物型偏微分方程

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 543
精华内容 217
关键字:

隐式求解偏微分方程