精华内容
下载资源
问答
  • Matlab偏微分方程快速上手:使用pdetool工具箱求解二维偏微分方程,适用于数学建模、数学实验,简单的偏微分方程数值计算与工程问题

    注:本人使用MatlabR2020a版本。

    1.pdetoolbox的调用

    打开MatlabR2020a,在命令行键入pdetool,进入pdetoolbox。
    输入pdetool进入pdetoolbox

    2.绘制定解区域(解的定义域)

    由图形界面可知,解的定义域是x,yx,y二维坐标构成的平面空间。我们必须设置自己的定解区域,才能定义自己的方程:
    导航栏下方的前5个按钮,分别对应绘制矩形求解区域、绘制按中心生成的矩形求解区域、绘制椭圆形(圆形)定解区域、绘制按中心生成的椭圆形(圆形)定解区域、绘制多边形求解区域。使用时,只需要点击后在绘图区域拖拽(多边形除外,多边形区域是在绘图区域点点以确定顶点),就可以生成定解区域了。
    这是前五个按钮
    上面这是前五个按钮。
    在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域。操作的时候用鼠标拖动操作柄拖拽就可以。(也可以画好几个叠起来)
    在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域
    真的是“随便”画一个就可以哦,因为在matlab下,求解区域的位置坐标精度达到了101610^{-16}左右,手动画几乎不可能画准。所以下一步教大家怎么细致地调节边界的坐标。

    3.手动调整定解区域的大小

    双击刚刚绘制好的区域,弹出一个对话框,里面是我们的定解区域的边界坐标信息(注意不全是坐标),我们可以在这里手动调整定解区域的位置(以矩形区域为例):
    刚刚打开时的样子
    这就是刚刚打开时的样子。因为这个区域是作者随便画的,所以坐标信息就像随机数一样。下面我们输入精确的数值:Left: -1, Bottom: -1, Width: 2, Height: 2, Name 就用默认的就好。
    输入参数

    参数的意义:Left:左边界的坐标(x左),Bottom:底边界的坐标(y底),Width:区域宽度,Height:区域高度。这样就得到了x[1,1],y[1,1]x\in[-1,1],y\in[-1,1]的矩形求解区域。调整好的求解区域显示效果如下。绘制好的区域

    4.调整绘图窗口的显示区域(调整显示坐标限)

    有时候我们会发现我们的定解区域太大了,绘图窗口显示不下;或者定解区域太小了,看上去非常不协调。上一个例子中作者的纵坐标显然非常吻合,但是横坐标多出来了(显示了左右两边的白框),那么我强烈建议大家调整完定解区域的坐标以后,再调整一下绘图窗口的显示区域。
    点击导航栏Options,再点击Axes Limits…(意为调整坐标限),可以手动设置坐标限。这里作者勾选了Auto,这样matlab将自动帮我们调整坐标限,使得定解区域位于界面中央。Options还有其他操作,大家可以自己尝试一下,这里就不介绍了。调整坐标限

    调整后的效果如下。

    调整效果

    5.确定边界条件

    点击导航栏下方第6个按钮(Ω\partial \Omega,意为Ω\Omega的边界条件),它用来显示边界。下面仍然以矩形区域为例。

    这个按钮
    双击边界

    此时显示了4个边界。双击任意一个边界,会弹出边界条件对话框,可以在这里随意设置边界条件。

    对话框

    在这里可以设置Dirichlet边界条件、Neumann边界条件和Robin边界条件(分别是第一、第二、第三边界条件),其中Robin边界条件和Neumann边界条件集成到一起了。按提示输入对应的系数就可。

    在这里作者使用了如下的边界条件:
    x=±1,u=0;y=±1,un=0.x=\pm1,u=0;y=\pm1,\frac{\partial u}{\partial n }=0.
    如果用了Dirichlet边界条件,边界将显示为红色;Neumann、Robin边界条件将显示蓝色,效果如下。

    边界条件

    6.确定偏微分方程的形式

    点击第7个图标(显示PDE字样),按提示输入偏微分方程的系数即可。在这里笔者求解波动方程:2u2t=u.\frac{\partial^2 u}{\partial^2 t}=\nabla u.

    在这里插入图片描述

    工具箱提供的方程通式如下:
    1.椭圆型Elliptic,通用数学形式为(cu)+au=f;-\nabla \cdot(c\nabla u)+au=f ;
    2.抛物型Parabolic,通用数学形式为dut(cu)+au=f;d\frac{\partial u}{\partial t}-\nabla \cdot(c\nabla u)+au=f ;
    3.双曲型Hyperbolic,通用数学形式为d2u2t(cu)+au=f;d\frac{\partial^2 u}{\partial^2 t}-\nabla \cdot(c\nabla u)+au=f ;
    4.特征值方程Eigenmodes,若λ\lambda为特征值,则数学形式为(cu)+au=λdu.-\nabla \cdot(c\nabla u)+au=\lambda d u .

    也可以自行指定求解的方程类型,比如比较常见的热传导、扩散等方程,可以在下面图示的下拉菜单中选择,但是仍然要按上面讲的方法手动设置系数。

    手动设置方程类型

    7.三角剖分

    由于Matlab pdetoolbox使用有限元方法求解,所以需要三角剖分。点击第8个图标(1个三角形图样)可以初始化剖分,点击第9个图标(4个三角形图样)可以增加剖分密度,这样可以提高计算精度,但是密度过高内存可能会爆掉,使用要谨慎。
    三角剖分

    8.设置初始条件,准备求解

    点击导航栏Solve,再点击Parameters…,进入求解参数设置器。在这里,第一行Time我们可以设置tt的求解范围及步长,默认情况下是不显示步长的(默认显示0:10意为从0求解到10,步长为1),我们按照Maltab等差数列的生成方法 a:j:b 就可以设置时间步长j了。

    第二行u(t0)、第三行u’(t0)表示t0t_0时刻的两个初始条件。这里作者使用了如下的初始条件。

    初始条件
    第四行和第五行表示相对容差和绝对容差,笔者查看了Matlab帮助中心,大概了解到这两个参数似乎与浮点数0的截断精度有关,太小的话会延长计算时间,如果你想了解更多,笔者把链接提供上来Absolute tolerance - MATLAB & Simulink - MathWorks 中国,假如我们对计算精度没有要求的话,使用默认值就可以了。这里笔者为了演示使用了0.001和0.0001。如果想跟着一起做,那么笔者把方程的代码也放上来:第一个是atan(cos(pi/2*x)),第二个是3*sin(pi*x).*exp(cos(pi*y))

    9.求解

    点击导航栏下方按钮(一个“==”字样的按钮,就是增加三角剖分密度右边那个按钮),这个按钮表示开始求解。如果求解完成的话会显示这个图。在这里可以点击“放大镜”按钮寻找感兴趣的区域放大来观察细节(放大之后想要缩小就要用上面步骤4的方法重新设置坐标限了,没有找到缩小的快捷键)。求解结果
    它直接显示了t=10t=10uuΩ求解区域\Omega的图像。这样的输出缺乏直观性,我们点击导航栏下方一个长得像matlab的logo的按钮(就是“==”按钮右边那个),调整绘图格式。

    输出格式窗口
    这个窗口有许多功能,作者就不再一一详述了。大家可以自行调试。比较常用的有“Contour”绘制等高线图,“Arrows”绘制向量场,“Height(3-D plot)”按3D模式输出(这个比较常用),“Animation”按动画形式输出(2D\3D都支持),此选项勾选后右边的Option选项会变亮,我们可以点进去在里面设置1秒显示的帧数、重复播放的次数。这里作者按照25阶变色、‘jet’ Colormap、向量场为u-\nabla u的形式静态输出t=0.5t=0.5时的结果,如下图所示。
    在这里插入图片描述

    这就是matlab pdetool工具箱的主要使用方法,本人也是小白一枚,所以欢迎大家批评指正,可以在评论区留下你的想法。

    参考:偏微分方程(姜礼尚《数学物理方程讲义》第三版)(更新完毕,附课件)来自于西北大学数统学院的马老师的数理方程视频课,是按数学系的讲法讲的数理方程,里面有那么两三个视频是讲如何用Matlab求解偏微分方程,如果你懂偏微分方程的话进去听一遍就会了。
    感觉应该是疫情期间这位老师的网课视频?那么西北大学的学生们也太幸福了,因为马老师讲的真的很好!人也很好玩哈哈,还去b站里别的老师的微分几何课程下面评论,正巧那个被评论的老师的同学就在b站讲拓扑哈哈,那个老师讲的也特别细我还给听完了(浙江理工庄老师)。扯远了!但是还是安利!!!

    展开全文
  • 先求通解再确定特解,是求常微分方程定解问题采用的方法,都某些偏微分方程,也能通过积分求出通解,进而确定出满足定解条件的特解。 两个自变量的一阶线性偏微分方程 今有两个自变量的一阶线性偏微分方程。 a(x,y)...

    先求通解再确定特解,是求常微分方程定解问题采用的方法,都某些偏微分方程,也能通过积分求出通解,进而确定出满足定解条件的特解。

    两个自变量的一阶线性偏微分方程

    今有两个自变量的一阶线性偏微分方程。
    a(x,y)ux+b(x,y)uy+c(x,y)u=f(x,y)(1) a(x,y)\frac{\partial u}{\partial x}+b(x,y)\frac{\partial u}{\partial y}+c(x,y)u=f(x,y) \tag{1}
    其中,系数a(x,y),b(x,y),c(x,y)a(x,y),b(x,y),c(x,y)是平面区域DR2D\subset \bold R^2上的连续函数,且a(x,y),b(x,y)a(x,y),b(x,y)不同时为零,f(x,y)f(x,y)在D上连续,称为方程的非齐次项。若f(x,y)0f(x,y)\equiv 0,方程为齐次的。

    思路:将两个自变量的方程化为求一个自变量的方程

    情况1:如果在D上,a(x,y)0,b(x,y)0a(x,y)\equiv 0,b(x,y)\neq 0,方程(1)改写为
    ux+c(x,y)b(x,y)u=f(x,y)b(x,y)(2) \frac{\partial u}{\partial x}+\frac{c(x,y)}{b(x,y)}u=\frac{f(x,y)}{b(x,y)} \tag{2}
    利用一阶线性常微分方程的求解方法得其通解。
    u(x,y)=exp(c(x,y)b(x,y)dy)[exp(c(x,y)b(x,y))f(x,y)b(x,y)dy+g(x)] u(x,y)=exp(-\int \frac{c(x,y)}{b(x,y)}dy)·[\int exp(\int \frac{c(x,y)}{b(x,y)})\frac{f(x,y)}{b(x,y)}dy+g(x)]
    其中,g(x)g(x)为任意C函数。

    情况2:如果在D上,a(x,y)b(x,y)0a(x,y)b(x,y)\neq 0,方程(1)不能直接积分求解,试作待定的自变量代换
    {ξ=φ(x,y)η=ψ(x,y)(a) \begin{cases} \xi=\varphi(x,y) \\ \eta=\psi(x,y) \end{cases} \tag{a}
    要求其雅可比(Jacobi)行列式
    J(φ,ψ)=(φ,ψ)(x,y)=φxφyψxψy0 J(\varphi,\psi)=\frac{\partial (\varphi,\psi)}{\partial(x,y)}= \begin{vmatrix} \frac{\partial \varphi}{\partial x} & \frac{\partial \varphi}{\partial y} \\ \frac{\partial \psi}{\partial x} & \frac{\partial \psi}{\partial y} \end{vmatrix} \neq 0
    以保证新变量ξ,η\xi,\eta的相互独立性,利用链式法则
    ux=uξξx+uηηx=φxuξ+ψxuηuy=uξξy+uηηy=φyuξ+ψyuη \frac{\partial u}{\partial x}=\frac{\partial u}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial u}{\partial \eta}\frac{\partial \eta}{\partial x}=\frac{\partial \varphi}{\partial x}\frac{\partial u}{\partial \xi}+\frac{\partial \psi}{\partial x}\frac{\partial u}{\partial \eta} \\ \frac{\partial u}{\partial y}=\frac{\partial u}{\partial \xi}\frac{\partial \xi}{\partial y}+\frac{\partial u}{\partial \eta}\frac{\partial \eta}{\partial y}=\frac{\partial \varphi}{\partial y}\frac{\partial u}{\partial \xi}+\frac{\partial \psi}{\partial y}\frac{\partial u}{\partial \eta}
    u=u(x,y)u=u(x,y)的方程(1)变为u=u(x(ξ,η),y(ξ,η))=u(ξ,η)u=u(x(\xi,\eta),y(\xi,\eta))=u(\xi,\eta)的新方程
    (aφx+bφy)uξ+(aψx+bψy)uη+cu=f(3) (a\frac{\partial \varphi}{\partial x}+b\frac{\partial \varphi}{\partial y})\frac{\partial u}{\partial \xi}+(a\frac{\partial \psi}{\partial x}+b\frac{\partial \psi}{\partial y})\frac{\partial u}{\partial \eta}+cu=f \tag{3}
    若取ξ=φ(x,y)\xi=\varphi(x,y)是齐次一阶线性偏微分方程
    a(x,y)φx+b(x,y)φy=0(4) a(x,y)\frac{\partial \varphi}{\partial x}+b(x,y)\frac{\partial \varphi}{\partial y}=0 \tag{4}
    的解,则新方程(3)称为(2)型的方程
    (aψx+bψy)uη+cu=f(b) (a\frac{\partial \psi}{\partial x}+b\frac{\partial \psi}{\partial y})\frac{\partial u}{\partial \eta}+cu=f \tag{b}
    η\eta积分便可求出通解。

    以下求解一阶线性偏微分方程(4),它的解对应与相应的常微分方程
    a(x,y)dyb(x,y)dx=0(5) a(x,y)dy-b(x,y)dx=0 \tag{5}
    亦即
    dxa(x,y)=dyb(x,y)(6) \frac{dx}{a(x,y)}=\frac{dy}{b(x,y)} \tag{6}
    的解之间存在确定的关系。

    定理:若φ(x,y)=h\varphi(x,y)=h(常数)是一阶常微分方程(5)在区域D内的隐式通解(积分曲线族),则ξ=φ(x,t)\xi=\varphi(x,t)是一阶线性偏微分方程(4)在区域D上的一个解。

    证明:设φ(x,y)=h\varphi(x,y)=h是方程(5)在D内的隐式通解,则过D内一点(x0,y0)(x_0,y_0)有一条积分曲线Γ0:φ(x,y)=φ(x0,y0)=h0\Gamma_0:\varphi(x,y)=\varphi(x_0,y_0)=h_0,此隐式解满足方程
    dxa(x,y)=dyb(x,y) \frac{dx}{a(x,y)}=\frac{dy}{b(x,y)}
    又沿此积分曲线Γ0\Gamma_0,有
    φxdx+φydy=0 \frac{\partial \varphi}{\partial x}dx+\frac{\partial \varphi}{\partial y}dy=0
    故在Γ0\Gamma_0上,有
    a(x,y)φx+b(x,y)φy=0 a(x,y)\frac{\partial \varphi}{\partial x}+b(x,y)\frac{\partial \varphi}{\partial y}=0
    由于(x0,y0)(x_0,y_0)是D内任意一点,故ξ=φ(x,y)\xi=\varphi(x,y)是一阶线性偏微分方程(4)在D上的解。

    定题的逆命题也成立。

    由常微分方程理论,一阶常微分方程(5)在区域D内存在且仅存在一族独立的积分曲线。如果求出了方程(5)的积分曲线族φ(x,y)=h\varphi(x,y)=h,再任取函数ψ(x,y)\psi(x,y),使在D上J(φ,ψ)0J(\varphi,\psi)\neq 0,以此φ\varphiψ\psi作变量代换(a)式,一阶线性偏微分(1)便可化为可积分求通解的方程(b)。

    特别地,当c(x,y)=f(x,y)0c(x,y)=f(x,y)\equiv 0时,方程(1)即为方程(4),相应的新方程(b)为uη=0\frac{\partial u}{\partial \eta}=0,其通解为u=g(ξ),g(ξ)u=g(\xi),g(\xi)为任意C函数。代回原自变量,得方程
    a(x,y)ux+b(x,y)uy=0 a(x,y)\frac{\partial u}{\partial x}+b(x,y)\frac{\partial u}{\partial y}=0
    的通解u=g(φ(x,y))u=g(\varphi(x,y))。这里,φ(x,y)=h\varphi(x,y)=h是常微分方程(5)的隐式通解,g(ξ)g(\xi)是任意C函数。

    如果给定u在某一曲线Γ:γ(x,y)=d\Gamma:\gamma(x,y)=d上的值,则需求解定解问题
    {a(x,y)ux+b(x,y)uy=0uγ(x,y)=d=θ(y) \begin{cases} a(x,y)\frac{\partial u}{\partial x}+b(x,y)\frac{\partial u}{\partial y}=0 \\ u|_{\gamma(x,y)=d}=\theta(y) \end{cases}
    用定解条件定出通解中的任意函数g(ξ)g(\xi)即可。

    这种求解定解问题的方法称之为通解法。常微分方程(5)或等价的方程(6)称为一阶线性偏微分方程(1)的特征方程,其积分曲线称之为特征曲线

    例1:求解右行单波方程的初值问题
    {ut+aux=0,t>0,<x<+ut=0=φ(x) \begin{cases} \frac{\partial u}{\partial t}+a \frac{\partial u}{\partial x}=0,\quad t>0,-\infty<x<+\infty \\ u|_{t=0}=\varphi(x) \end{cases}
    其中,a>0a>0为常数。

    :特征方程dxadt=0dx-adt=0特征线族xat=hx-at=h

    ξ=xat,η=x\xi=x-at, \eta=x,则方程化为
    uη=0 \frac{\partial u}{\partial \eta}=0
    η\eta积分得通解u=g(ξ)=g(xat)u=g(\xi)=g(x-at),其中,g(ξ)g(\xi)是任意C函数。由初始条件
    ut=0=g(x)=φ(x) u|_{t=0}=g(x)=\varphi(x)
    得该初值问题的解u(t,x)=φ(xat)u(t,x)=\varphi(x-at)

    (x,u)(x,u)平面上看(图1.3.1),t=0t=0u=φ(x)u=\varphi(x),对每个固定时刻t>0u=φ(xat)]t>0,u=\varphi(x-at)],其图形相当于曲线u=φ(x)u=\varphi(x)向右移动了atat,波形的传播速度为aa。称这样的解为右行波解

    (x,t,u)(x,t,u)空间看(图1.3.2),波形沿特征线传播。在特征线xat=hx-at=h上,u=φ(xat)=φ(h)u=\varphi(x-at)=\varphi(h),故当观察者沿某条特征线前行时,看到的波形始终不变。如果初始扰动φ(x)\varphi(x)只发生在区间h1xh2h_1\leq x\leq h_2内,则这个扰动沿着(x,t)(x,t)平面上的特征条形域h1xath2h_1\leq x-at\leq h_2传播。

    在这里插入图片描述

    本例中初始条件给在非特征线的直线t=0t=0上,从通解可唯一确定特解。如果初始条件给在一条特征线xat=h0x-at=h_0上,初始条件uxat=h0=g(xat)xat=h0=g(h0)=φ(x)u|_{x-at=h_0}=g(x-at)|_{x-at=h_0}=g(h_0)=\varphi(x),当φ(x)\varphi(x)为非常值函数时无解,当φ(x)\varphi(x)为常数φ0\varphi_0时有无穷多个解g(ξ)g(\xi),只要g(h0)=φ0g(h_0)=\varphi_0

    对于左行单波方程utaux=0\frac{\partial u}{\partial t}-a\frac{\partial u}{\partial x}=0,同样可求得其通解为左行波u=g(x+at)u=g(x+at)gg为任意C1C^1函数(一阶连续导数)。

    展开全文
  • 偏微分方程(PDE) 偏微分方程理论: 物理/工程问题————翻译(建模)/物理工程规律————》数学问题(PDE) 物理/工程问题————求解/数学理论————》数学结果 物理/工程问题————分析————》...

    1. 简介

    微分方程:描述自然界中存在的物理现象和普遍规律。

    • 常微分方程(ODE)
    • 偏微分方程(PDE)

    偏微分方程理论:

    • 物理/工程问题————翻译(建模)/物理工程规律————》数学问题(PDE)
    • 物理/工程问题————求解/数学理论————》数学结果
    • 物理/工程问题————分析————》数学公式/物理意义

    偏微分方程的基本概念:

    • 定义:未知函数及其偏导数所满足的方程;F(x, u(x), Du, D2u,…, Dnu) = 0;
    • 阶数:偏微分方程中偏导数的最高阶数,有n阶就为n。
    • 线性偏微分方程:方程中的未知函数或其偏导数的项都是一次项
    • 齐次方程:每项都含有未知函数或未知函数(定义的u = u(x, y))的偏导数。

    导出微分方程的基本步骤:

    1. 明确研究的物理量
    2. 明确物理量遵守的物理规律
    3. 按物理规律写出微分方程(泛定方程)
    4. 微元法:划出一个微元,分析临近部分和ta的相互作用

    微分方程解决的主要问题:

    • 描述对象特征随时间变化的过程
    • 分析对象特征的变化规律
    • 预报对象特征的未来形态
    • 研究控制对象特征的手段

    常用的物理定律概述

    1. 牛顿第二定律:F = ma;
    2. 胡克定律,在弹性限度内,弹性体的张应力和弹性体的形变量成正比:张应力 = 杨氏模量 * 相对伸长,
    3. 热传导的傅立叶定律:在dt时间内,通过面积元dS流入一个小体积元的热量dQ,与沿着面积元法线方向的温度梯度成正比,也与dS和dt成正比:(k是导热系数,由材料决定)
    4. 牛顿冷却定律:单位时间内从周围介质,传到边界上单位面积的热量,与表面和外界的温度差成正比:(这里u1是外界媒质的温度,h为比例系数)
    5. 扩散定律(斐克定律):单位时间流过某横截面的杂质量dm与该横截面积S和浓度梯度成正比:(式中D为扩散系数,负号是指杂质浓度在减少)

    2. 栗子——香烟过滤嘴的作用

    问题:

    • 过滤嘴的作用和ta的材料与长度有什么关系;
    • 人体吸入的毒品量和哪些因素有关,因素的影响程度怎么样

    模型分析:

    • 分析吸烟时毒物进入人体的模型,建立吸烟过程模型;
    • 假设一个机器人持续吸烟,并且环境不变

    模型假设(假设的变量都在这儿哦):

    • 烟草长 l1,过滤嘴长 l2,总长度 l = l1 + l2,假设毒品均匀分布;
    • 毒物进入空气和沿香烟穿行的数量比:a’ : a, a’ + a = 1;
    • 未点燃烟草和过滤嘴的吸收率(过滤程度)分别为b和ß;
    • 烟雾沿着香烟的穿行速度为常速v,香烟的燃烧速度为u,v >> u;
    • Q为毒物进入人体的总量
    • q(x, t)为毒物的流量,w(x, t)为毒物密度

    模型建立:

    • 建立毒物进入身体的总量:Q = ƒ0T q(l, t)dt, T = l1 / u;

    • 求q(x, 0) = q(x)——流量守恒

      q(x) - q(x + ∆x) = bq(x)∆∂, 0 ≤ x ≤ l1;

      q(x) - q(x + ∆x) = ßq(x)∆∂, l1 ≤ x ≤ l2;

      (∆∂ = ∆x / v)

      同时假设:

      • q(0) = aH0;(香烟毒物被吸入的量)
      • H0 = uw0;(香烟毒物的减少量)

      可得到

      dq/dx = -b/v * q(x), 0 ≤ x ≤ l1;

      dq/dx = -ß/v * q(x), l1 ≤ x ≤ l2;

      &

      q(x) = aH0e-bx/v, 0 ≤ x ≤ l1;

      q(x) = aH0e-bx/ve-ß(x-l1)/v, l1 ≤ x ≤ l2;

    • 目标是求q(l, t),假设t时刻,香烟燃烧到了x = ut;

      将q(x)中的H0拓展为H(t),x随时间会变成x-ut,l1会变成l1 - ut

      可以求得:

      q(x, t) = aH(t)e-b(x-ut)/v,ut ≤ x ≤ l1;

      q(x, t) = aH(t)e-b(l1-ut)/ve-ß(x-l1)/v,l1 ≤ x ≤ l;

    • 求w(ut, t)求∆t内毒物密度的增量。

      w(x, t + ∆t) - w(x, t) = b(q(x, t))/v * ∆t(单位长度烟雾被吸收的部分)

      我们已知:

      • q(x, t) = aH(t)e-b(x-ut)/v;
      • H(t) = uw(ut, t);(H0的拓展公式)

      结果可得到:

      w(ut, t) = w0/a’ * (1 - ae-a’but/v), a’ = 1 - a;

    • 计算Q,吸入一致烟毒物进入人体的总量;

      根据求出来的 w(x, t) 和 q(l, t) 可以代入Q = ƒ0T q(l, t)dt, T = l1 / u;求出现在的Q = aM e-ßl2/v µ®, r = a’bl1/v, µ® = 1 - e-r / r;

    结果分析:

    1. Q和a,M成正比,aM是毒物集中在x = l处的吸入量;

    2. e-ßl2/v是过滤嘴的因素;

    3. µ®是烟草的吸收作用;

    4. 判断与不带过滤嘴的香烟比较

      Q1和Q2的差别为b和ß;

      我们已知:

      • 带滤嘴:Q1 = (aw0v / a’b) e-ßl2/v(1 - e-a’bl1/v);
      • 不带滤嘴:Q~12 = (aw0v / a’b) e-bl2/v(1 - e-a’bl1/v);

      所以明显是滤嘴提高了吸收能力

    特点:

    • 引入两个基本函数:流量q(x,t)和密度w(x,t),运用物理学的守恒定律建立微分方程,构造动态模型。

    3. 偏微方程的导出

    3.1. 波动方程的导出

    • 传输线方程(电报方程 )
    • 均匀薄膜的微小横振动方程
    • 流体力学与声学方程
    • 电磁波方程

    ta们的物理本质根本不同,但这些方程的数学形式与弦振动方程和杆纵振动方程完全一样:

    3.2. 扩散方程的导出

    3.2.1. 细杆热传导

    现象描述:导热细杆各点的温度是不均匀的,热量由高到低传导。

    目的:求出细杆中温度的变化方程。

    定律:

    • 热传导的傅立叶定律:q = -k∆u(单位时间内流过单位时间的热量q与温度的梯度成正比,k为热传导系数);
    • 热量守恒定律:热量变化 = 边界流入 + 热源放出;

    第一种情况(系统无热源):

    (热传导仅由物体内部温度不均引发)

    u(x, t)为x处t时刻的温度。

    一维情况下热传导的傅立叶定律有: q = -k∂u/∂x;(q为热流强度)

    在x方向上(微元法):

    • dt时间流入左表面的热量为:q|x Adt;
    • dt时间流出右表面的热量为:q|x+dx Adt;
    • 所以净流入为:q|x Adt - q|x+dx Adt = -∂q/∂xAdxdt =(然后把q换掉再整合)= kuxxAdxdt;

    因为热量Q = c(øAdx)[u|t - u|t+dt] = c(øAdx)utdt;(假设ø为密度)

    所以更具热量守恒定律得到kuxxAdxdt = c(øAdx)utdt;

    结果:ut - a2uxx = 0, a2 = k/cp(这就得到了均匀物体内部无热源时的热传导方程)

    第二种情况(系统内有热源):

    比第一问不过就是多个热源,也就是加一个热量呗(设F(x, t)为单位时间,单位体积内产生的热量:kuxxAdxdt + F(x, t)Adxdt = c(øAdx)utdt;

    得到:ut - a2uxx = F(x, t) / cp = f(x, t);(ƒ(x, t)为热源强度)

    三维情况下就是:ut - a2∆u = f(x, t);

    3.2.2. 扩散问题

    扩散现象:系统浓度不均用时,物质从高浓度转移向低浓度的现象。

    目的:建立空间各点浓度u(x, y, z, t)的方程。

    物理规律:扩散定律和粒子数守恒定律

    • 扩散定律(裴克定律):单位时间内流过单位面积的分子数或质量,与浓度的梯度成正比:(D为扩散系数)

      沿着n方向的大小:

    • 粒子数守恒定律:单位时间流入一定体积粒子数流出同一体积粒子数的差,等于该体积单位时间内粒子数的增加量。

      同样是微元法,划出一个小立方体v分析浓度变化规律。

      • 计算单位时间沿x-方向的净流入量(负号是表示和浓度梯度的方向相反):(利用两个平面的流入流出的积分)
      • 同理可以求出y方向和z方向的净流入量:
      • 所以总的净流入量为
      • 单位时间的粒子增加数为:
      • 根据粒子数守恒定律就可以得到:

      我们在前面得到过这样的结果:

      所以我们将扩散定律代入就可以得到三维扩散方程

      我们484看到了吗有很多D,如果扩散是均匀的,D就是一个常数,我们可以令D = a2,则有

      那如果这个空间里面有扩散源,也就是扩散是从这儿来的,那么我们可以将ta们相等得到公式:

    3.3. 稳定场方程

    3.3.1. 热传导方程

    如果边界条件与热源内不随时间变化,一定时间后,内部温度会达到稳态。

    u = u(x, y, z), ∂u/∂t = 0;

    so

    ut - a2∆u = f(x, y, z, t) -> ∆u = -f(x, y, z)泊松方程

    3.3.2. 静电场下的泊松方程和Laplace方程

    从高斯定理出发,我们可以结合扩散定律简单的得到:

    我们可以得到

    还因为

    我们可以得到两个方程:

    • 泊松方程:
    • Laplace方程(当的时候):

    4. PDE导出总结

    系统的物理状态除了取决于自己状态,还取决于系统环境的状态。

    • 物理、工程在t>0时刻的系统环境,在数学中称为边界条件;
    • 物理、工程在t=0时刻的系统状态,在数学中称为初始条件;

    定解条件:是边界条件和初始条件的总体,反映了问题的个性。

    泛定方程:不带边界和初始条件的方程,反映问题的个性。

    4.1. 初始条件 ——描述系统的初始状态

    • 波动方程:含有时间的二阶导数,所以需要两个初始条件
      1. u|t=0 = u(x) 系统各点的初始位移
      2. ∂u/∂t|t=0 = v(x) 系统各点的初始速度
    • 热传导方程:含有时间的一阶导数
      1. u(x, t)|t=0 = u(x)初始时刻的温度分布
    • 泊松方程和拉普拉斯方程:不含时间的倒数,不需要初始条件

    注意:初始条件给的是初始状态下物理量的分布,而不是指某一位置

    4.2. 边界条件

    未知函数在边界满足条件:

    1. 第一类Dirchlet边界条件:已知未知函数在边界上的函数值:

      基本形式:u(x, y, z, t)| = ƒ(M, t);

    2. 第二类Neumann:已知未知函数在边界上法线方向的导数值;
      基本形式:∂u(x, t)/∂n|x0 = ƒ(x0, t);

    3. 第三类Robin混合边界条件:混合牛顿冷却定律、傅立叶实验定律(一维);

      牛顿冷却定律:

      q = h(u - ø)n,u为物体表面温度,ø是周围介质温度,h是热交换系数。一维条件下简化为q = h(u - ø)。

      牛顿冷却定律可以得出流出热量与外界温差成正比。


      傅立叶实验定律:

      q|x=a = -k∂u/∂x|x=a

      基本形式:(u + H∂u/∂n)| = ƒ(x0, y0, z0, t);

    4.3. 小结

    5. 偏微分方程的求解

    5.1. 达朗贝尔公式:行波解

    基本思想:

    1. 求PDE的通解;
    2. 用定解求特解

    关键步骤:通过变量变换,将波动方程化为齐次二阶偏微分方程。(所以也叫做波动方程的初始问题,或者柯西问题)

    求解问题:

    我们根据式1可以拆分得到,变换变量也就可以转换为

    通过复合求导我们可以得到:
    ——》

    一些不要的东西删了:

    求通解:

    对两个创建的变量进行积分:

    积分得到

    再积分得到

    我们的通解就是

    其中ƒ1, ƒ2 是连续可微的一元函数。

    ƒ1(x - at)和ƒ2(x + at)的意义?

    • x - at为正方向运动的行波
    • x + at为反方向运动的行波

    确定待定函数形式,求特解:

    我们已经有初始条件:中的后两条。

    我们可以根据之前的条件可以进行两个公式的化简:

    1. image-20210714150847740
    2. image-20210714150921073

    所以可以得到我们求出的达朗贝尔公式:
    image-20210714151011495

    5.2. 分离变量法

    基本思想:将偏微分方程通过分离变量变成一个常微分方程!

    关键步骤:求特征值问题

    适用问题:有界域上的波动方程、热传导方程、稳定场方程等

    求解问题:

    特征值问题:含有未知常数的常微分方程,求非零解的问题;

    特征值:是方程有非零解的常数值;

    特征函数:和特征值相对于的非零解

    5.3. Fourier变换法(傅立叶变换法)

    基本思想:将偏微分方程通过Fourier变换变成一个常微分方程~

    关键步骤:求常微分方程定解问题和ta的解的方法,Fourier逆变换

    适用问题:无界域上的波动方程、热传导方程等

    求解问题:

    Fourier变换法基本步骤:

    1. 对偏微分方程与初始条件中实行傅立叶变换,转化为常微分方程;
    2. 解常微分方程的定解问题,得到相应的傅立叶变换式;
    3. 对该式进行逆傅立叶变换,求的原问题的解
    展开全文
  • 偏微分方程的数值解(一):定解问题 & 差分解法 偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法 偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布 偏微分...

    偏微分方程的数值解系列博文:

    偏微分方程的数值解(一):定解问题 & 差分解法

    偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

    偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

    偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

    偏微分方程的数值解(五): 二维状态空间的偏微分方程的 MATLAB 解法

    偏微分方程的数值解(六): 偏微分方程的 pdetool 解法


    目录

    1 方程类型

    (i)椭圆型偏微分方程                                 (ii)抛物型偏微分方程          

    (iii)双曲型偏微分方程                                (iv)特征值问题

    (v)非线性椭圆偏微分方程                         (vi)方程组

    2 边界条件

    (i)Dirichlet 条件:​                     (ii)Neumann 条件: ​      

    (iii)对于偏微分方程组,混合边界条件为

    3 求解偏微分方程

    例 6 求解泊松方程                                                                 例7 考虑最小表面问题

    例8 求解正方形区域上的热传导方程                                  例9 求解正方形区域上的波方程

    例 10  求解泊松方程                                                           习题


    MATLAB 中的偏微分方程(PDE)工具箱是用有限元法寻求典型偏微分方程的数 值近似解,该工具箱求解偏微分方程具体步骤与用有限元方法求解偏微分方程的过程是 一致的,包括几个步骤,即几何描述、边界条件描述、偏微分方程类型选择、有限元划 分计算网格、初始化条件输入,最后给出偏微分方程的数值解(包括画图)。 下面我们讨论的方程是定义在平面上的有界区域Ω 上,区域的边界记作∂Ω 。

    1 方程类型

    MATLAB 工具箱可以解决下列类型的偏微分方程:

    (i)椭圆型偏微分方程

    其中c, a, f 和未知的u 可以是Ω 上的复值函数。

    (ii)抛物型偏微分方程

    其中c,a, f ,d 可以依赖于时间t 。

    (iii)双曲型偏微分方程

    (iv)特征值问题

    其中λ 是未知的特征值,d 是Ω 上的复值函数。

    (v)非线性椭圆偏微分方程

    其中c, a, f 可以是u 的函数。

    (vi)方程组

    2 边界条件

    边界条件有如下三种:

    (i)Dirichlet 条件:

    (ii)Neumann 条件:

    这里 \overrightarrow{n} 为区域的单位外法线,h,r, q, g 是定义在∂Ω 上的复值函数。

    对于二维方程组情形,Dirichlet 边界条件为    

    Neumann 边界条件为:

    (iii)对于偏微分方程组,混合边界条件为

    这里 μ 的计算是使得满足 Dirichlet 边界条件。

    3 求解偏微分方程

    例 6 求解泊松方程 

                          求解区域为单位圆盘,边界条件为在圆盘边界上u = 0。

    解 它的精确解为

    下面求它的数值解,编写程序如下:

    %(1)问题定义
    g='circleg'; %单位圆
    b='circleb1'; %边界上为零条件
    c=1;a=0;f=1;
    %(2)产生初始的三角形网格
    [p,e,t]=initmesh(g);
    %(3)迭代直至得到误差允许范围内的合格解
    error=[]; err=1;
    while err > 0.01,
    [p,e,t]=refinemesh(g,p,e,t);
    u=assempde(b,p,e,t,c,a,f); %求得数值解
    exact=(1-p(1,:).^2-p(2,:).^2)/4;
    err=norm(u-exact',inf);
    error=[error err];
    end
    %结果显示
    subplot(2,2,1),pdemesh(p,e,t);
    subplot(2,2,2),pdesurf(p,t,u)
    subplot(2,2,3),pdesurf(p,t,u-exact')

    例7 考虑最小表面问题

    g='circleg';
    b='circleb2';
    c='1./sqrt(1+ux.^2+uy.^2)';
    rtol=1e-3;
    [p,e,t]=initmesh(g);
    [p,e,t]=refinemesh(g,p,e,t);
    u=pdenonlin(b,p,e,t,c,0,0,'Tol',rtol);
    pdesurf(p,t,u) 

    例8 求解正方形区域上的热传导方程

    求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的热传导方程

    边界条件为Dirichlet条件u = 0。

    解 这里是抛物型方程,其中c = 1, a = 0, f = 0, d = 1。编写程序如下:

    %(1)问题定义
    g='squareg'; %定义正方形区域
    b='squareb1'; %边界上为零条件
    c=1;a=0;f=0;d=1;
    %(2)产生初始的三角形网格
    [p,e,t]=initmesh(g);
    %(3)定义初始条件
    u0=zeros(size(p,2),1);
    ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);
    u0(ix)=1
    %(4)在时间段为0到0.1的20个点上求解
    nframe=20;
    tlist=linspace(0,0.1,nframe);
    u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);
    %(5)动画图示结果
    for j=1:nframe
     pdesurf(p,t,u1(:,j));
     mv(j)=getframe;
    end
    movie(mv,10) 

    例9 求解正方形区域上的波方程

                  求解正方形区域{(x, y) | −1 ≤ x, y ≤ 1}上的波方程

    解 这里是双曲型方程,其中 c = 1, a = 0, f = 0, d = 1。编写程序如下:

    %(1)问题定义
    g='squareg'; %定义正方形区域
    b='squareb3'; %定义边界
    c=1;a=0;f=0;d=1;
    %(2)产生初始的三角形网格
    [p,e,t]=initmesh(g);
    %(3)定义初始条件
    x=p(1,:)';y=p(2,:)';
    u0=atan(cos(pi*x));
    ut0=3*sin(pi*x).*exp(cos(pi*y));
    %(4)在时间段为0到5的31个点上求解
    n=31;
    tlist=linspace(0,5,n);
    uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);
    %(5)动画图示结果
    for j=1:n
     pdesurf(p,t,uu(:,j));
     mv(j)=getframe;
    end
    movie(mv,10) 

    例 10  求解泊松方程

    求解区域为单位圆盘,边界条件为在圆盘边界上u = 0。

    解 它的精确解为

    下面求它的数值解,编写程序如下:

    g='circleg';
    b='circleb1';
    c=1;a=0;f='circlef';
    [p,e,t]=initmesh(g);
    [p,e,t]=refinemesh(g,p,e,t);
    u=assempde(b,p,e,t,c,a,f);
    exact=-1/(2*pi)*log(sqrt(p(1,:).^2+p(2,:).^2));
    subplot(2,2,1),pdemesh(p,e,t);
    subplot(2,2,2),pdesurf(p,t,u)
    subplot(2,2,3),pdesurf(p,t,u-exact') 

    习题

    1. 求二维拉普拉斯方程

    2. 求初边值问题

    在0 ≤ t ≤ 3范围内的数值解。

    3. 求热传导方程初边值问题


    偏微分方程的数值解系列博文:

    偏微分方程的数值解(一):定解问题 & 差分解法

    偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

    偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

    偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

    偏微分方程的数值解(五): 二维状态空间的偏微分方程的 MATLAB 解法

    偏微分方程的数值解(六): 偏微分方程的 pdetool 解法

    展开全文
  • n个自变量的一阶线性偏微分方程(n≥2n\geq2n≥2) n个自变量的一阶线性偏微分方程的一般形式为 ∑j=1nbj∂u∂xj+cu=f(7) \sum_{j=1}^nb_j\frac{\partial u}{\partial x_j}+cu=f \tag{7} j=1∑n​bj​∂xj​∂u​+cu...
  • 偏微分方程-特征线法(转)

    万次阅读 2012-06-03 10:55:57
    只要初始不是沿着特征线给定,即可通过特征线法获得偏微分方程的精确解。 其基本思想是通过把双曲线型的准线性偏微分方程转化为两组常微分方程,再对常微分方程进行求解。两组常微分方程中的一组用于定义特征线...
  • 偏微分方程基础

    2021-05-25 09:31:28
    第50篇 偏微分方程 一个偏微分方程(PDE)包含两个或多个自变量的导数。这与之前描述的常微分方程(ODE)相反,后者只涉及一个自变量。 工程和科学中的许多现象都是用偏微分方程描述的。例如,一个因变量,如压力或温度...
  • 在数列、微分方程、矩阵三个不同的领域中都见到了“特征值”这个名字,难道仅仅是因为他们都刻画了问题的特征吗?当时我想,一定不是的。或许这些特征值是同一个特征值。因为本人水平较低,这个...
  • MATLAB偏微分方程数值解视频课程

    千次阅读 2019-05-17 10:07:44
    结合MATLAB偏微分方程数值解工具箱介绍偏微分方程的求解,分GUI和MATLAB函数两种实现方式进行介绍。 【课程收益】 MATLAB偏微分方程数值解工具箱的使用 有限单元法 用GUI和MATLAB编程两种方式求解PDE问题 第一...
  • 图 像处理作为一种预处理的手段,几乎成为所有图像处理方法的前奏。在许多情况下,图像滤波...满足对比度不变和仿射不变的偏微分方程只有一个,即 AMSS(Affine Morphological Scale Space)方程.L.Alvarez,F.Guichar...
  • 来源:新浪了凡春秋的博客在科学技术各领域中,有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程偏微分方程,再加上边界条件、初始条件构成的数学...
  • 高等数学 偏微分方程 波动方程
  • 偏微分方程
  • COMSOL 偏微分方程接口

    千次阅读 2020-06-05 19:03:58
    更多阅读:sppy.site COMSOL Multiphysics是一款强大的多物理场仿真软件。除了内置的众多典型物理场外,其还允许用户基于方程进行建模,例如使用 PDEs (Partial Differential...COMSOL的偏微分方程接口下常用的有两种.
  • 偏微分方程(Partial Differential Equation I) 偏微分方程(Partial Differential Equation II) 格林函数法 积分变换法 非线性数学物理问题 参考文献: 《数学物理方法》| 吴崇试 《数学物理方法》| 梁昆淼 《数学...
  • 图片来源-维基百科小角度简单摆若最高处( )的绳子和最低处(速度最大)的绳子的角度为 ,则可使用下列公式算出它的振动周期。 公式证明摆球受力分析绳与对称线夹角为 ,绳长为 ,绳距离对称线的水平距离为 。对于...
  • 通过非线性的偏微分方程扩大梯度空间、保留梯度较大的边缘,增强图像的纹理细节。由于多光谱图像阴暗波段的纹理较弱,不容易辨别其所有信息,为了更好地使增强效果完全体现出来,使用直方图均衡化来调节亮度的不...
  • 二阶线性偏微分方程

    万次阅读 2016-08-22 23:01:46
    这三类方程的形状很特殊,它们是二阶线性偏微分方程的三个典型代表。 一般形式的二阶线性偏微分方程之间的共性和差异,往往可以从对这三类方程的研究中得到。 一维弦振动方程是双曲型的,一维热传导方程是抛物...
  • 针对具有微小缺失或破损的图像修复问题,基于变分和偏微分方程理论基础提出一种TV双调和偏微分方程图像修复方法。该模型在Sobolev对偶分布空间H-1(Ω)上考虑一个变分问题,它的Euler-Lagrange方程是一个与多孔介质...
  • 应用偏微分方程 作 者: 谷超豪,李大潜,沈玮熙 著 出版时间:2014 丛编项: 高等学校教材 内容简介  本书的写作意图是通过几个经过选择的主题的简单介绍,使读者了解偏微分方程应用的一些基本内容和特点,以增强...
  • 本书系统讲解偏微分方程及其定解问题的求解方法,通过大量实例讨论偏微分方程解的性质。 目录 第8章 拉普拉斯变换和汉克尔变换及其应用 第9章 有限差分数值方法 第10章 抽样和离散傅里叶分析及其在偏微分方程中的...
  • 和涉及微分方程特征值问题。 与传统的 ODE 求解方法相比,在目标函数足够平滑的情况下,谱方法自然具有收敛速度超快的优势。 有关光谱方法的更多详细信息,请查看 。 它列出了用于理解谱方法和 MATLAB 项目Chebfun...
  • 这些讲座中的 ipython 笔记本涵盖了使用局部和全局搭配求解偏微分方程的基本方面。 该系列方法包括标准有限差分、使用 Chebyshev、Hermite、Laguerre、傅立叶和 Sinc 基函数的谱方法,以及具有各向同性误差的多维...
  • 某些二阶线性偏微分方程,可分解为两个一阶线性偏微分方程,有可能积分求出通解。例如,二阶方程 ∂2u∂x∂y+∂u∂x=0 \frac{\partial^2u}{\partial x\partial y}+\frac{\partial u}{\partial x}=0 ∂x∂y∂2u​+∂...
  • 椭圆型偏微分方程数值解法

    千次阅读 多人点赞 2020-08-04 19:33:11
    % 一维椭圆方程求解(常微分方程值问题) % -u'' + q(x)u = f(x), 0<x<1, 取q(x) = x, f(x) = (x-1)exp(x) % u(0) = 1, u(1) = e; 边界条件 % 真解为 u = exp(x) N = 20; h = 1/N; x_all = (0:h:1)'; x = x_...
  • 实质上是一种信号的滤波器,其用途是信号的平滑处理,我们知道数字图像用于后期应用,其噪声是的问题,由于误差会累计传递等原因,其实编程运算的话就是一个模板运算,拿图像的八连通区域来说,中间点的像素就等于...
  • 前言 本文是作者对本学期所学《数学物理方程》和《数值分析》中涉及偏微分方程求解的知识所进行的一个类型归纳和方法总结,并配合一定数目的例题进行相关练习。 ...
  • >利用Matlab求解二阶偏微分方程的一般有以下步骤    → 题目定义:由方程(3.4.33)和(3.4.35)可以看出,参量 是二阶偏微 分方程的主要参量,只要这几个参量确定,就可以定下偏微分方程的结构。此外要做的事是确定偏...
  • 从通解求特解的方法对多数偏微分方程定解问题不适用。求解线性偏微分方程定解问题较多采用以下做法:先求方程的一族特殊解,再将这族特殊解作适当线性叠加求得定解问题的特解。接下将介绍这种做法的理论根据——线性...
  • 近十年多来,偏微分方程(PDEs)的理论和方法在图像处理各领域的应用越来 越引起了人们的关注。本文在现有基于偏微分去噪模型的基础上,对模型的线性扩散、 非线性扩散和保真项进行了具体的研究,结合“自蛇”(self...

空空如也

空空如也

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

偏微分方程特征值问题